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Introduction générale 


INTRODUCTION GÉNÉRALE 


Les véhicules aériens autonomes sont, de nos jours, des moyens de plus en plus utilisés. Ils sont 
aussi nommés drones, vecteurs où UAVs. Le fait que le véhicule soit autonome (ou commandé à 
distance) permet d'effectuer des missions longues (la fatigue du pilote n’est pas un élément d’arrêt 
de mission) et en toute sécurité pour le pilote. 


Ces véhicules autonomes sont classés en plusieurs catégories, en fonction de leurs missions res- 

pectives et de leurs possibilités de charges utiles (l'équipement pouvant être embarqué) [97, 123] : 

e Les drones tactiques (TUAV). Ce sont des drones plutôt petits. Ils n’ont pas une portée d’action 
élevée et sont assez peu encombrants. L’autonomie de ce type de drone varie de quelques 
minutes à quelques heures. 

e Les drones volant à moyenne altitude et de longue endurance (MALE). Ces drones permettent 
l’utilisation d’une charge utile pouvant atteindre 500 kg. Leur altitude de croisière est comprise 
entre 5 000 et 15 000 mètres. L’autonomie de ce type de drone peut être comprise entre 20 et 
40 heures. 

e Les drones volant à haute altitude et de longue endurance (HALE). L’altitude de croisière de 
ces drones est supérieure à 15 000 mètres. 

Les drones de combat (UCAV) sont des engins pouvant embarquer une charge utile létale. Ils 

peuvent être tactique, MALE ou HALE. 


Dans ce travail, nous ne nous intéresserons qu'aux drones HALE et MALE. De plus nous consi- 
dérerons que la charge utile n’est autre qu’une caméra. En particulier, nous nous intéresserons aux 
problèmes de planification de trajectoires pour le système composé du vecteur et de sa charge utile. 


0.1 Contexte 


Le travail présenté dans ce document est le fruit d’une thèse financée par un FUI (Fonds Unique 
Interministériel), dirigée par le Professeur J.-P. Gauthier (Université du Sud-Toulon-Var) et co- 
dirigée par J.-F. Balmat (Université du Sud-Toulon-Var) et U. Serres (Université Claude Bernard 
Lyon 1). 


Les drones constituent une technologie de plus en plus mise en avant dans différents projets 
[133, 77]. Le laboratoire LSIS se joint à cette recherche grâce à cette thèse subventionnée par le pro- 
jet SHARE. Supporté par un consortium d'entreprises !, le projet a pour objectif la création d’une 
station de contrôle au sol universelle pour des drones à voilure fixe ou tournante (e.g., des avions ou 
des hélicoptères). 


Le but du laboratoire LSIS est d’assurer le couplage entre la trajectoire du drone (le vecteur) et 
ses capteurs (la charge utile). Ce couplage vecteur/capteur définit l’asservissement de la trajectoire 
du vecteur à la manipulation de la charge utile par l'opérateur, i.e. nous souhaitons conserver le 
champ de vision de l’opérateur en optimisant un des paramètres suivants : 

e la discrétion, 


1. Opéra Ergonomie, ONERA, Thales Alénia Space, Eurocopter, Adetel group 
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e l'angle de dépression de la caméra pour disposer de la meilleure observation (—45 à —30 
degrés), 

e l'orientation du soleil (éviter l'orientation face au soleil, éviter l'ombre portée du vecteur), 

e la maximisation des capacités de zoom. 


Dans certaines situations, le couplage vecteur/capteur asservira la manipulation de la caméra 
(et non celle du vecteur). C’est le cas, par exemple, lors de l’observation d’une zone prédéfinie via 
un modèle donné pour la caméra (e.g., nous souhaitons faire un état des lieux de tout un terrain). 


Afin de satisfaire au mieux à cet objectif, nous devons prendre en compte les contraintes sui- 
vantes : 
e Contraintes environnementales : 
hauteur du sol, 
sens du vent, 
obstacles, 
proximité du relief, 
couloirs aériens imposés (Flight levels), 
couloirs aériens interdits (No fly zone), 


O 


ontraintes sur le vecteur : 

domaine de vol, 

paramètres de vol (voilure fixe ou tournante), 

paramètres du vecteur (ce que l’utilisateur peut commander), 


ontraintes sur le capteur : 
libertés de rotation, 
vitesses de rotation, 
temps de rafraîchissement, 


oo 0 © O0 0 6 6° 45 0 6. © OC 6 0 


En plus des contraintes citées précédemment, le vecteur et les capteurs sont soumis à des modes 
de pilotage différents (choisis par l'utilisateur) et agissant à des niveaux de contrôle distincts [95, 94] : 
e Modes de commandes du vecteur : 
o Mode de base : l'utilisateur commande le cap, la vitesse de vol, la vitesse d’ascension. 
o Mode automatique : l’utilisateur demande au vecteur de naviguer vers un point de passage 
défini. 
o Mode pattern : l'utilisateur demande au vecteur de suivre un pattern [91], ï.e. une figure 
prédéfinie comme un cercle, un 8, un hippodrome.….. 
e Modes de commandes du capteur : 
o Mode manuel : l’utilisateur commande le gisement, le site, le niveau de zoom. 
o Mode automatique 1 : l'utilisateur peut choisir un asservissement sur un point fixe. 
o Mode automatique 2 : l'utilisateur peut choisir un asservissement sur une cible mouvante. 


INTRODUCTION GÉNÉRALE 0.2. Contribution 


Ainsi, toutes ces questions d’asservissement du vecteur et de la charge utile conduisent à une 
problématique de planification de trajectoires contraintes par les éléments cités précédemment. Cette 
planification de trajectoires doit permettre d’effectuer les missions demandées. Elles peuvent être : 

e Observer un point fixe en utilisant un pattern autour du point d'intérêt à altitude constante, 

e.g., un bâtiment dont nous souhaitons connaître les caractéristiques, 

e Observer un point fixe sans utiliser de pattern autour du point d’intérêt, e.g., une maison, une 

entrée d’une grotte ou d’une cache, 

e Observer un point mobile dont la trajectoire est connue, e.g., un convoi ami à protéger, 

e Observer un point mobile dont la trajectoire est inconnue, e.g., un convoi de chars ou de 

fantassins, 

e Effectuer un balayage d’une zone d'intérêt fixe avec la caméra, 

e Effectuer un balayage d’une surface mobile, i.e. effectuer une observation libre autour d’un 

point mobile, e.g., une protection de convoi. 


Afin de satisfaire au mieux à ces missions, les algorithmes de planification devront permettre au 
vecteur et à la caméra de : 
e Observer un point d'intérêt à 360 degrés, 
e Observer un point d'intérêt sous un secteur angulaire (relatif à la cible), 
e Observer un point d'intérêt sous un angle (azimut et élévation) souhaité. Ceci implique, entre 
autre, la capacité de spécifier l’angle d'observation par l’opérateur via l'interface utilisateur 


(IHM). 


À terme, grâce au projet SHARE, l'utilisateur pourra effectuer simplement toutes sortes de 
missions, Comme : 

e Surveiller un objectif, 

e Renseigner sur un objectif, 

e Rechercher un objectif, 

e Acquérir un objectif, 

e Évaluer les dommages sur un objectif, 

e Effectuer un relais de communication. 


Remarque 0.1. Pour certaines missions, l’objectif peut être soit ami, soit ennemi. 


Le but de ce travail est donc de construire des algorithmes de planification de trajectoires pour 
le vecteur permettant de prendre en compte la charge utile ainsi que les contraintes inhérentes 
au problème. Les algorithmes présentés par la suite sont voués à être utilisés en amont du mode 
automatique de commande, i.e. les algorithmes développés fournissent des points de passage à suivre 
par le vecteur et sa charge utile. 


0.2 Contribution 


Dans la suite, nous supposerons que la charge utile est une caméra positionnée sous le vecteur. 
De plus, nous considérons que le vecteur est soit un drone de type HALE (Haute Altitude Longue 
Endurance) volant à vitesse et altitude constante, soit un drone MALE (Moyenne Altitude Longue 
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Endurance) volant à vitesse constante. 


Nous verrons que les équations du mouvement du drone HALE sont alors celles données par le 
système de la voiture de Dubins ? [51], i.e. 


à = cos 0 
ÿ = sin 0 
= u, 


avec (x,y,0) € R? x S1 l’état du système (où (x, y) € R? sont les coordonnées du drone dans le plan 
d'altitude constant, et 9 le cap), et u est la variable de contrôle. Ce système a beaucoup été étudié, 
notamment dans [92, 116, 117]. 


Pour un drone de type MALE, nous verrons que les équations du mouvement peuvent être 
considérées comme celles de l’évolution de son repère cinématique [141, 71], i.e. 
X = RV 
R= RQ, 
avec (X,R) € R° x SO(3) l’état du système où X représente les coordonnées du drone dans le 
repère global et À est la matrice de rotation entre le repère cinématique du drone et le repère global. 
Vo = (1,0,0) représente le vecteur vitesse unitaire dans le repère du drone et est dans l’alignement 
du corps de l’appareil. La matrice Q représente la vitesse angulaire du solide considéré dans le repère 


cinématique du drone. C’est une matrice antisymétrique contenant les variables de contrôle u1, u2 
et u3 correspondant, respectivement, aux vitesses angulaires du roulis, tangage et lacet du drone. 


Pour les modèles présentés ci-dessus, nous allons, dans un premier temps, proposer des méthodes 
de planification afin de répondre aux contraintes décrites au paragraphe précédent. 

Dans un deuxième temps, nous répondrons à une problématique qui est naturellement apparue 
durant ce projet : comment proposer des trajectoires comparables à celles qu’un pilote expérimenté 
serait amené à décrire s’il était directement aux commandes de l’appareil? L'objectif étant de se 
conformer au mieux aux règles de l'Organisation de l'Aviation Civile Internationale [74, circulaire 
328, article 2.13]. 


2. Le système de Dubins est utilisé pour modéliser simplement une voiture. Nous l’utilisons pour notre vecteur en 
considérant une altitude fixe. 
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0.3 Présentation du manuscrit 


Le manuscrit est organisé comme suit : 
e La partie [| propose des méthodes de planification et répond au problème de couplage vec- 
teur/capteur : 


o Le chapitre 1 introduit le problème de planification. 


o L'état de l’art, la présentation des modèles choisis ainsi que des premières réponses au pro- 
blème de planification sont proposées en chapitre 2. 


o Le chapitre 3 est une présentation de notre contribution avec les principaux résultats. 
o Le chapitre 4 dresse le bilan de la partie I. 

e La partie II traite la problématique anthropomorphique. Dans cette partie, nous proposons 
une méthode, basée sur le contrôle optimal inverse, permettant au drone d’agir comme le ferait 
un pilote : 


o Nous introduisons cette problématique au chapitre 5. 


o Le chapitre 6 présente un état de l’art des méthodes bio-inspirées et énonce le problème de 
contrôle optimal inverse. 


o Les résultats obtenus, pour un drone HALE, sont présentés dans le chapitre 7. 


o Le chapitre 8 rappel les points importants de la partie Il. 


Première partie 


Planification de trajectoires 


Chapitre 1 


Introduction 


Au sein du projet SHARE, notre travail est d’assurer le couplage entre la trajectoire du drone (le 
vecteur) et ses capteurs (la charge utile). Ceci passe notamment par de la planification de trajectoires. 

Dans cette partie, nous verrons quelles sont les méthodes actuellement utilisées pour répondre à 
ce type de problématique. Nous présenterons aussi celles que nous avons utilisées. 


1.1 Le problème de planification 


Le problème de planification soulevé a deux facettes. D'abord nous voulons que le drone se 
déplace selon une trajectoire d’un point à un autre de l’espace. Ensuite, nous voulons qu’en fonction 
de la trajectoire du drone, la charge utile embarquée remplisse sa mission. 

Ce problème, aussi appelé couplage vecteur/capteur, implique l’utilisation de méthode de pla- 
nification de trajectoires pour obtenir un chemin permettant de réaliser simultanément ces deux 
objectifs en ne compromettant pas les différentes contraintes inhérentes aux deux parties. 


Le but est donc de trouver un (ou plusieurs) algorithme de planification pour que le vecteur et 
la charge utile atteignent leurs objectifs respectifs. 

Par exemple, nous demandons au vecteur de se diriger vers une position fixe (en évitant les 
obstacles, ou les zones aériennes interdites) tout en demandant à la caméra embarquée de suivre une 
cible (sans la perdre de vue). 

Le problème considéré peut être posé de la manière suivante : 


Le problème de couplage vecteur/capteur 


Étant donnée une position initiale du vecteur et de sa caméra (notée qo), comment calculer une 
trajectoire q(t), du drone et de sa caméra, définie sur [0,7] (où T est un paramètre fixé ou libre) 
qui permette d’atteindre un point voulu &cible- 


La réponse à cette question sera donnée par la suite. Nous nous intéressons à deux cas de figure : 
un point final fixe (e.g., rejoindre un point de passage ou un pattern prédéfini) ou mobile (e.g., suivre 
un convoi). 
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Dans ces deux cas, il faut calculer une trajectoire de référence notée q(t), définie sur [0,T|, sa- 
tisfaisant les équations de la dynamique du vecteur et de la caméra, joignant les points qo et qcible 
spécifñés tels que q(0) = qo et q(T) = Gcible- 


La planification de trajectoires se fait en général en optimisant un critère : minimisation des 
efforts du moteur d’une voiture, de la consommation de carburant d’une fusée, maximisation du 
confort des passagers, minimisation de la longueur du chemin parcouru, minimisation du temps 
pour relier deux points... 

De telles trajectoires optimisant un critère peuvent être obtenues grâce au Principe du Maximum 
Pontryagin, noté PMP [3]. Une application de cette méthode est présentée dans le paragraphe 2.2.3. 

Dans de nombreuses applications, le PMP est utilisé en amont du développement pour trouver 
des morceaux de trajectoire (aussi appelés primitives) optimaux localement. La trajectoire sera alors 
simplement un recollement de primitives [70, 72]. 


Une autre méthode, très utilisée, consiste à discrétiser l’espace [57, 137, 96] et à appliquer des 
algorithmes de cheminement (basés sur la théorie des graphes), tels que les algorithmes de Dikjstra 
[1, 136], A* [66, 140] ou d’autres encore [50]. 


D’autres manières de procéder existent, pour une présentation plus complète, voir [83]. 


Pour répondre au mieux à certains types de missions, la norme en vigueur [95] définit des trajec- 
toires types, appelées patterns. Par exemple, lorsque le drone doit protéger un convoi, il doit suivre 
une trajectoire adaptée qui lui permet de surveiller le convoi et son voisinage proche. 

Ces patterns sont aussi définis et utilisés dans les rapports [91, 72]. Les patterns souvent utilisés 
pour aider au bon déroulement des missions sont : 

e Le cercle : ce pattern permet de faire de la surveillance d'objectif, 

e L’hippodrome : ce pattern permet de protéger un convoi, 

e Le huit : ce pattern permet de bombarder un objectif. 


Lors de l’utilisation de patterns pour répondre aux besoins d’une mission, la question principale 
est de déterminer comment le vecteur peut le rejoindre. Après une étude des méthodes couramment 
utilisées, au paragraphe 2.3, nous en proposerons de plus élaborées au paragraphe 3.2. 


Chaque méthode de planification présentée est développée dans le but d’être utilisée dans un 
module haut niveau, voué à être exécuté sur une machine distante. Celui-ci doit permettre de trans- 
mettre des points de passage à un module d’auto-pilotage du drone, voir Figure 1.1 ou, pour plus 
de détails, [4, Figure 5] (cf. Annexe B), [6, Figure 4] (cf. Annexe A), [25, Figure 3]. 

Cette architecture est définie par la norme STANAG 4586 [95] et elle est résumée dans l’article 
[16]. La méthode décrite dans [100] permet de créer un module bas niveau sur le drone, i.e. un 
auto-pilote. Dans [108], les auteurs utilisent un module de pilotage automatique du commerce. 


Dans le cas d’une cible mobile, i.e. pour la poursuite de convoi, cette architecture peut être 
modifiée pour fournir toutes les données nécessaires à l’algorithme de planification. 
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Module de . 
pilotage automatique 


rl Actionneurs || Capteurs 
a Régulateur | 


Station au sol | 


Interface = Module de 


utilisateur! planification 
) 


FIGURE 1.1 — Schéma de l’architecture considérée 


En effet, si la cible est mobile, alors le point final gp est une fonction du temps. Cette fonction 
du temps doit être connue à chaque instant pour que l'algorithme puisse calculer une trajectoire. 
Quand la cible mobile est amie, cette information est supposée connue à chaque instant. Lorsque la 
cible est ennemie, ces données ne sont pas toujours connues. 

Pour pallier ce manque d’information, nous devons réaliser une estimation de la trajectoire. Le 
schéma de l’architecture correspondant à ce cas de figure est représenté dans [25, Figure 3]. 


Rappelons que le problème de couplage vecteur/capteur doit aussi prendre en compte les contraintes 
énoncées en Introduction Générale. Nous devons mettre en place des méthodes permettant de ne 
pas entrer en conflit avec celles-ci. 


1.2 Organisation de la partie I 


Dans cette partie, nous traitons du problème de planification de trajectoires. Le chapitre 2 
concerne l’état de l’art de la planification de trajectoires. Tout d’abord nous considérons le cas où 
la cible est un point fixe et dans un second temps la cible sera mobile. 


Ensuite, dans le chapitre 3, nous présentons nos principaux résultats. Dans ce chapitre, le para- 
graphe 3.2 est le plus important : il décrit nos méthodes de planification vers un pattern donné. 

L'étude du couplage vecteur/capteur est faite au paragraphe 3.1. Cette étude montre que le 
problème de couplage peut se restreindre à l’étude de la planification du vecteur seul. 

Enfin, le paragraphe 3.4 propose une application des méthodes proposées pour effectuer du suivi 
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de convoi en milieu encombré. 
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CHAPITRE 2. ÉTAT DE L'ART 2.1. La modélisation 


Dans ce chapitre, nous présentons quelques méthodes de planification de trajectoires pour un 
drone, dans le cas du suivi d’une cible fixe ou mobile. 

Tout d’abord, nous abordons la modélisation des drones. La plupart des méthodes sont basées 
sur des modèles de comportement cinématiques des drones. 

Pour répondre au problème de planification point-point, i.e. les points initiaux et finaux sont 
fixés, nous donnons ensuite un aperçu de méthodes de stabilisation et de contrôle. 

Pour terminer, nous proposons des généralisations de ces méthodes pour effectuer du suivi de 
cible mobile. La cible est considérée comme amie ou ennemie, i.e. sa trajectoire future est connue 
ou non, et évolue dans un milieu encombré ou non. 

Nous ne présentons ici que les méthodes de planification usuelles. Pour une bibliographie plus 
complète, nous vous renvoyons aux livres de LaValle [83] et Laumond [82]. Un aperçu des méthodes 
de suivi de convoi est présenté dans les articles [143, 42]. 


2.1 La modélisation 


Nous allons présenter certains modèles cinématiques. [ls sont fréquemment utilisés en planifica- 
tion de trajectoires de véhicules évoluant dans un plan ou dans l’espace usuel à 3 dimensions. 


Dans la suite nous considérons des systèmes de la forme générale 
d= fau), qEeQ, ue% (21) 


où q est l’état du système et Q une variété lisse de dimension n. Le contrôle est donné par u € Y 
et f(-,u) est un champ de vecteur lisse, i.e. C®, sur la variété Q. 


Nous présentons les équations de chaque modèle, ainsi que les propriétés de contrôlabilité asso- 
ciées. 

Dans la plupart des cas, la contrôlabilité est intuitivement évidente. Pour rappel, elle est définie 
comme suit [15] : 


Définition 2.1 (Contrôlabilité d’un système). Si, pour tout couple (qo; Gb) de deux états du 
système (2.1), il existe un temps T et une fonction localement intégrable u : [0,T] — 7 telle que la 
solution du problème 


ÿ = f(q,u), q(0) — 40 


satisfasse q(T) = qible, alors le système de contrôle est dit contrôlable. 


Cette notion de contrôlabilité d’un système est une propriété importante. Par exemple, si le 
système est non contrôlable, la question d’existence de solutions optimales se pose. 


Jean-Michel Coron, dans son cours [41], présente quelques méthodes d’analyse de la contrôla- 
bilité. Il cite, pour les systèmes linéaires, le critère de Kalman qui donne une condition nécessaire 
et suffisante de contrôlabilité d’un système. Par contre, dans le cas de systèmes non-linéaires, la 
contrôlabilité n’est pas toujours évidente à prouver. 
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FIGURE 2.1 — Schéma d’un drone HALE 


Pour les systèmes non-linéaires et sans drift, le théorème de Rashevski-Chow ([15, Théorème 2.3] 
ou [126]) donne une condition nécessaire de contrôlabilité partout locale, i.e. en temps petit il existe 
un contrôle qui permet d’atteindre tout point d’un certains voisinage du point de départ, sans en 
sortir. 


2.1.1 Modélisation des drones HALE 


Nous présentons une modélisation des drones de type HALE (Haute Altitude Longue Endu- 
rance). Puisque ces appareils volent à vitesse et altitude constantes, nous considérons l’étude de 
leurs mouvements dans leurs plans d'altitude fixe. 

Un schéma d’une situation générique du drone HALE dans le plan est donné en Figure 2.1. Sur 
ce schéma, (x,y) définissent les coordonnées du centre de gravité du drone dans le plan d’altitude 
constante et 0 définit son cap (la flèche représente le vecteur vitesse du système). 

Un modèle largement utilisé pour ce type de véhicules est le système de Dubins [51, 16] 


à = cos 0 
ÿ = sing (22) 
Ô = u, 


en notant q = (x,y,0) € R? x S1 l’état du système et u la variable de contrôle (généralement 
contrainte). 


Ce système est de la forme générale (2.3), dite affine dans le contrôle : 
4 = f{g) + ug(q), (2.3) 
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avec 
x cos Ô 0 
q= {y}, f(a)={|snm0|, g(a)= 10 
0 0 1 


Remarque 2.2. La commande u appliquée au modèle de Dubins (2.2) représente la courbure de la 
trajectoire suivie par le système. En effet, rappelons que la courbure & d’une courbe paramétrée 
dans un repère orthonormé y(t) = (x(t),y(t)) s'exprime par 
_ y 
E (à? + ÿ2)3/2° 
Le long d’une courbe du système (2.2), nous avons 
_ Ô cos? 0 + 6 sin? 0 
(cos? 0 + sin? 9)?/2 
= U. 
Dans le cas des drones HALE, la commande « est contrainte par le rayon de courbure minimal du 
vecteur. L'ordre de grandeur de ce rayon de courbure est le kilomètre. 


Remarque 2.3. Dans la modélisation (2.2), nous considérons que les drones HALE ont une vitesse 
normalisée égale à 1. Parfois (voir le paragraphe 3.2) nous considérons uw € [—-1,1]. Ces deux hypo- 
thèses sont possibles grâce à une normalisation du modèle cinématique du drone HALE. 

En effet, si nous considérons que le drone HALE a une vitesse constante v et une commande 
vérifiant u € [—Uuimax; Umax], alors le système cinématique peut s’écrire sous la forme 


& = vcos 0 
ÿ = vsin 0 (2.4) 
ô=u. 
En utilisant la normalisation en temps et espace suivante 
rl 
Ut 
TL =UT 
Ÿ = wy, 


avec La = L/Umax Et Lb — Umax/V, le système normalisé vérifie (2.2) avec un contrôle compris dans 
l'intervalle [—1, 1]. 


Un système similaire à (2.2) a été étudié par J.A. Reeds et L.A. Shepp [110]. Il autorise le véhi- 
cule à faire marche arrière, i.e. vu € {—1,1} dans le système (2.4). 


La contrôlabilité des systèmes de Dubins et de Reeds-Shepp est en fait évidente. La raison 
théorique pour laquelle ces systèmes sont contrôlables est exposée dans le livre de Steven M. LaValle 
[83]. Pour une étude complète, se référer à [127]. 
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2.1.2 Modélisation des drones MALE 


Le paragraphe 3.2 traite de la problématique de la planification dans le cas d’un drone MALE 
(Moyenne Altitude Longue Endurance) qui évolue dans l’espace à 3 dimensions. 


Dans ce cas, le modèle considéré doit représenter les principales caractéristiques d’un drone 
MALE, i.e. nous souhaitons pouvoir modifier sont cap ainsi que son altitude. Nous supposons de 
nouveau qu’il a une vitesse constante unitaire. 

Considérons un système de Dubins (similaire au drone HALE) que nous étendons en associant 
une commande indépendante sur l’altitude, comme suit 


à = cos 0 
ÿ = sin 0 

2.5 
_. (2.5) 
Ô = up, 


avec q = (x,y,z,0) € R° x S!, l'état du système (où (x,y,z) € R° sont les coordonnées du drone 
dans l’espace, et 9 son cap). Les variables de contrôle sont alors u9, pour commander le cap, et u,, 
pour modifier l'altitude. 

Le modèle ainsi formé est appelé modèle de type avion de Dubins et a été utilisé de nombreuses 
fois dans certains problèmes de planification [10, 94, 60]. 


Ce modèle ne représente pas l’évolution du drone en tant que corps solide dans R?. 

Certains modèles de drones MALE d’un niveau de complexité plus élevé sont présentés dans 
[26, 58, 62]. 

Si nous restons au niveau purement cinématique, les équations classiques de la mécanique du 
solide constituent une représentation raisonnable d’un drone MALE. Ces équations constituent 
d’ailleurs la base des modèles plus précis [26, 58, 62]. 


Nous considérons donc l’évolution du repère de l’avion (le repère cinématique) par rapport au 
repère global (galiléen, nous supposons la terre plate et immobile). 

L'état q associé au drone est composé de la matrice de passage du repère cinématique au repère 
global ainsi que des coordonnées X = (x1,%2,x3) du centre de masse dans le repère global (cf. Figure 
D). 

En reprenant les notations de la Figure 2.2, l’état peut s’écrire sous la forme q = (X, R) avec la 
matrice de rotation À vérifiant 


ed — 


Li J F] = R Fe Juav Kuav 
Nous supposons que le repère cinématique est contrôlé en vitesse angulaire par trois commandes 


correspondant aux commandes du roulis, tangage et lacet du drone. Pour plus de détails, se référer 
à [26, p.43]. 
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Au final, le système cinématique peut s’écrire sous la forme [141, 71] 
. 
2 (2.6) 
R=RQ 
avec 
eq;—=(X,R)E R° x SO(3), l'état du système où X représente les coordonnées du drone dans le 
repère global et R est la matrice de rotation entre le repère cinématique du drone et le repère 
global, 
e V5 — (1,0,0), le vecteur vitesse normalisé à 1 dans le repère du drone, qui se trouve dans 
l’alignement du corps de l’appareil, 


O0 —u3 u2 ; x à TT à 2 
e ( — us 0 ) € 50(3), la vitesse angulaire du solide considérée dans le repère cinéma- 
—U2 U1 


tique du drone (50(3) est l’algèbre de Lie de SO(3), constituée des matrices antisymétriques). 


Les variables de contrôle u1, u2 et u3 correspondent, respectivement, aux vitesses angulaires du 
roulis, tangage et lacet. 


La contrôlabilité du modèle (2.6) a été étudiée d’un point de vue mathématique dans [119] et est 
obtenue en appliquant [115, Théorème 6.7] (ou bien [28, Théorème 1]). Elle est encore intuitivement 
évidente. 


Remarque 2.4. En ne considérant que des rotations autour de l’axe k et des commandes de roulis et 
tangage nulles, l’évolution des trajectoires du système (2.6) n’est rien d’autre que celle du système 
de Dubins (2.2) dans un plan horizontal d'altitude fixée dès que la vitesse initiale de l’avion est dans 
un plan horizontal. 


2.2 Planification point-point 


Disposant maintenant de représentations mathématiques des véhicules considérés, nous sommes 
en mesure de traiter le problème de planification. Nous allons, dans cette partie, étudier quelques 
méthodes applicables dans le cas d’une cible fixe. Nous appelons ce problème la planification de 
trajectoires point-point. 

De nombreuses méthodes sont utilisées pour répondre à ce problème. Citons, par exemple, les 
travaux de Jean et al. [75, 38, 88] basés sur la théorie du contrôle ou celui de Takeiï et al. [128] qui est 
basé sur une résolution numérique des équations de Hamilton-Jacobi. Pour une bibliographie plus 
étendue, se référer à [83]. 


2.2.1 La platitude 


Une première méthode utilisable pour répondre à ce problème de planification de trajectoires est 
l’utilisation de la propriété de platitude des systèmes. 

La notion de platitude a été développée par M. Fliess et al. au début des années 1990 [55]. Elle 
est présentée dans un article des Techniques de l'ingénieur [111], inspiré du livre [90], de la manière 
suivante : 
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NS 


To 


“T1 


= _, 


O | , 


FIGURE 2.2 — Schéma du drone en 3D 


“Cette propriété, qui permet de paramétrer de façon très simple le comportement dy- 
namique d’un système, est basée sur la mise en évidence d’un ensemble de variables 
fondamentales du système : ses sorties plates. |... 

Une sortie plate est composée d’un ensemble de variables qui permet de paramétrer 
toutes les autres variables du système, l’état, la commande, mais également la sortie.” 


M. Fliess, dans [51|, définit la platitude comme suit 
Définition 2.5 (Platitude d’un système). La platitude est une propriété d’un système d'équations 


différentielles ordinaires, sous-déterminé, c’est-à-dire à moins d'équations que d’inconnues. Le sys- 
tème est dit plat si il existe m variables y = (y1,...,Ym) telles que : 


1. Toute variable du système s'exprime comme fonction différentielle de y, i.e. une fonction des 
composantes de y et de leurs dérivées jusqu’à un ordre fini. 


2. Toute composante de y s'exprime comme fonction différentielle des variables du système. 
3. Les composantes de y et leurs dérivées sont linéairement indépendantes. 
y est appelé sortie plate ou linéarisante. 


L'intérêt de la platitude est qu’elle permet une forme de linéarisation exacte du système, ce qui 
permet de se ramener à un problème classique de commande linéaire. 


Dans le domaine de l’ingénierie, la commande par platitude est surtout utilisée pour les pro- 
blèmes de poursuite de trajectoires. M. Fliess et al. présentent de nombreux exemples d'applications 
dans les articles [55, 54, 53]. Cette méthode a été reprise pour commander l’entrée dans l’atmosphère 
d’une navette spatiale [93] et pour réaliser des fonctions de suivi de trajectoires en temps réel [138]. 
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Plus récemment, elle a été utilisée pour mimer la préhension des rapaces par un quadcopter [131]. 


Bien que cette méthode soit utilisée assez fréquemment, nous ne l’appliquerons pas pour faire 
de la planification de trajectoires. Par contre, nous pouvons, avec cette méthode de contrôle, mettre 
en œuvre une loi de guidage pour effectuer du suivi de trajectoire, i.e. s’en servir dans le module 
de pilotage automatique. Dans ce cas de figure, une méthode de planification doit être utilisée pour 
fournir des points de passage, comme proposé dans [131, Figure 7]. 


2.2.2 La méthode de Lyapunov : une utilisation du principe de LaSalle 


La méthode de Lyapunov est une méthode de stabilisation vers un point ou un ensemble de 
points choisis (e.g., un pattern). 


Avant d'expliquer cette méthode, rappelons la définition de la stabilité. 


Définition 2.6 (Stabilité au sens de Lyapunov [102]). En reprenant les notations précédemment 
définies, considérons un système différentiel de la forme q = f(q), q € Q. 

Supposons qu’il existe un point d'équilibre q € Q, i.e. f(q) = 0. 

L'équilibre g € Q est dit stable (au sens de Lyapunov) si et seulement si pour tout € > 0, il existe 
n > 0 tel que pour toute condition initiale qo vérifiant ||go — g|| < mn, la solution de 4 = f(q) issue de 
go à t — 0, est définie pour tout temps t positif et vérifie [[g(t) — q| < €. 

Si le système n’est pas stable, il est dit instable. 


En plus de la définition de la stabilité autour d’un point d'équilibre, nous définissons la stabilité 
asymptotique et la stabilité globalement asymptotique, pertinente dans notre cas. 


Définition 2.7 (Stabilité asymptotique [102]). En reprenant les notations de la Définition 2.6, 
l'équilibre q est dit localement asymptotiquement stable si et seulement si il est stable et si, de 
plus, il existe 7 > 0 tel que toutes les solutions q(t) de 4 — f(q), partant en t — 0 de conditions 
initiales qo telle que ||go — g|| < n, convergent vers q lorsque t tend vers +oo. 


Lorsque la condition initiale peut être librement choisie, g est dit globalement asymptotique- 
ment stable. 


Le théorème de Lyapunov ou le principe d’invariance de LaSalle [80], qui est une extension 
de celui-ci, sont appliqués pour prouver la stabilité globalement asymptotique d’un système sans 
contrôle q = f(q). Ils sont énoncés dans [102, 135] et d’une manière plus concise dans [61, Théorème 
3.5, p.190]. Nous les rappelons ci-après. 


Définition 2.8 (Fonction de Lyapunov [135]). Nous reprenons les notations de la Définition 2.6. 
Soit Q un ouvert de R” contenant le point d'équilibre g. La fonction V : Q — R est une fonction 
de Lyapunov en g sur Q si 

e V'est Cl sur (1, 

e V(g) =0 et pour tout ge Q\{q}, V(q) > 0, 
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e V décroît le long des trajectoires de g = f(q), i.e. pour tout q € Q, 


VG) = (VV f@) = 2,7 (DA < 0 


Si, pour tout qg € Q \ {q}, V(q) < 0, la fonction de Lyapunov est dite stricte. 


Théorème 2.9 (Théorème de Lyapunov [135]). Nous reprenons les notations de la Définition 2.8. 
Si il existe une fonction de Lyapunov au point d'équilibre q sur Q, alors le point q est stable. 
Si la fonction de Lyapunov est stricte alors q est asymptotiquement stable. 


Théorème 2.10 (Principe d’invariance de LaSalle [135]). Nous reprenons les notations de la défi- 
nition 2.6. 
Soit V : Q — R+ une fonction de classe C\ telle que 


e V est propre sur Q, i.e. VT E V(Q), V-{([0,T]) est compact dans Q, 
° VyEQ,  (VV(g), f(g)) < 0. 


Soit T le plus grand sous-ensemble de {q € Q | (VV (a), f(q)) = 0} invariant par le flot (en temps 
t>0) de q = f(q). Alors toute solution q(t) de q = f(q) tend vers T. 


Remarque 2.11. Dans le Théorème 2.10, si l'se réduit à un point d'équilibre q alors q est globalement 
asymptotiquement stable dans (. 


P. Rouchon et al. expliquent [102, Partie 3.5.2] comment utiliser ces résultats pour stabiliser le 
système (2.3). 

La méthode, dite de Lyapunov, consiste à synthétiser, en se servant du Théorème 2.9 ou du 
Théorème 2.10, une loi de commande qui conduit à la stabilité asymptotique globale. 


Cette méthode est utilisée dans de nombreuses applications. Par exemple, certains auteurs ap- 
pliquent cette méthode pour éviter les collisions entre les drones [17] ou bien pour commander des 
robots bio-inspirés [36]. 

Le problème de suivi de convoi sera aussi traité par cette méthode, cf. paragraphe 2.3.4. 


Nous nous inspirons de cette méthode au paragraphe 3.2 pour le problème de planification point- 
pattern (1.e. le drone doit rejoindre un ensemble de points). 


La méthode de Backstepping 


Une autre méthode de stabilisation, présentée dans [76, Paragraphe 14.3], est construite à partir 
du principe de LaSalle. Elle s’applique à une certaine classe de systèmes non-linéaires à un contrôle, 
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appelée strict-feedback form par l’auteur, et donnée par 


& = fo(x) + go(x)21 
à = f(x, A) + g(x, 1)22 


22 = abs 22) in g2(x, 21; 22)23 


da = fe Pier) Fo Mises 10) 
de Ce 


où æ € R” et 21 à zx sont des scalaires et u € R le contrôle. De plus, dans cette classe de sys- 
tèmes, nous supposons que les f;, à € {0,...,k}, s’annulent à l’origine et que g(x,21,...,2;) £ 0, 
i € {0,...,k} sur tout le domaine. 


En utilisant itérativement le principe de LaSalle, nous construisons une commande uw qui permet 
de stabiliser ces types de systèmes. Cette méthode de stabilisation est dite de backstepping. 


Une utilisation de la méthode de backstepping sur un système simple est proposée dans [102, 
Paragraphe 3.5.2]. Dans cet exemple, en reprenant les notations précédentes, l’auteur considère 
k=1et à = u. 


2.2.3 Le problème de contrôle optimal 


Les méthodes de contrôle optimal sont utilisées pour piloter les drones tout en minimisant une 
certaine fonction de coût (appelée critère), notée J. Ce critère à optimiser peut, par exemple, être 
le temps de parcours ou la longueur de la trajectoire. 

L'outil classique pour résoudre ce type de problème est le Principe du Maximum de Pontryagin, 
noté PMP [3]. Nombreuses sont les utilisations de ce principe. La mécanique hamiltonienne en est un 
exemple où le critère à optimiser n’est autre que l’action, i.e. l’intégrale sur le temps du Lagrangien. 
Dans ce cas, ce principe est dit d’action minimale et la quantité de mouvement fait office de ce que 
nous nommerons, par la suite, vecteur adjoint [2]. 


Dans la suite, nous noterons ( P;) le problème optimal associé à une fonction de coût L et nous 
le définissons de la manière suivante : 


Problème de contrôle optimal (P) 


Minimiser le coût intégral J(q,u) = 1 L(q(t),u(t))dt parmi tous les contrôles admissibles u(-) qui 
permettent, sous la contrainte du système (2.1), de relier un point initial go à un point final qcible 
en un temps 7’, libre ou fixé. 


Remarque 2.12. Un point qcibie étant donné, l’ensemble des extrémales solutions du problème (P;) 
à partir de tout point qo est appelé synthèse optimale du problème (P7). 
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Pour répondre au problème ( P;), il faut d’abord étudier l’existence des ses solutions. Générale- 
ment, l’utilisation du théorème de Fillipov [3, Théorème 10.1] permet de prouver un tel résultat. 

Ensuite, les trajectoires optimales sont obtenues en utilisant le Principe du Maximum de Pon- 
tryagin, dont voici l'énoncé [19] 


Théorème 2.13 (Principe du Maximum de Pontryagin (PMP)). Considérons le système de contrôle 
de la forme (2.1) et un coût de la forme J(q,u) = Fe L(q(t),u(t))dt où L est une fonction lisse par 
rapport à q. Le problème est de relier un point initial qo à un point final qeinie en un temps T fixé. Si le 
couple (q(-),u(-)) est optimal sur [0, T], alors il existe une application p(-) : [0, T] — R” absolument 
continue appelée vecteur adjoint, et un réel À < 0, tels que, pour tout t € [0, T] le couple (X,p(t)) est 
non trivial, et tels que, pour presque tout t € [0,T]|, 


(2.7) 


où A(X,p,q,u) = (p, f(q,u)) + AL(q,u) est le Hamiltonien du système, et de plus, nous avons la 
condition de maximisation presque partout sur [0,T] 


HA PE), q(E), u(t)) = max Æ(A,p(6), q(t),v) = Cste (2.8) 
vEY 
Remarque 2.14. e Si de plus le temps final T' pour rejoindre la cible qinie n’est pas fixé, nous 


avons la condition suivante : .Æ#7(X, p(t), q(t),u(t)) = 0, Vte [0,T]. 

e Il existe des versions plus générales du PMP {[3, 139]. 

e Ce théorème est une condition nécessaire d’optimalité. 

e Si le contrôle n’est pas borné ou qu'il se trouve à l’intérieur du domaine Z des contrôles 
admissibles, alors l'équation (2.8) entraîne DL (Xp, qu) = (. 

e Définissons une trajectoire extrémale comme un quadruplet (À,p(-),q(-),u(-)) qui vérifie les 
conditions nécessaires d’optimalité du PMP. 
Si À — 0, l’extrémale est dite anormale. 
Si À < 0, l’extrémale est dite régulière. Dans ce cas, le vecteur adjoint étant défini à un scalaire 
multiplicatif près, nous pouvons choisir les valeurs usuelles À = —1 ou À = —1/2. 

e Si les points q0, Gaible ne sont pas fixés mais appartiennent à des variétés lisses Qo et Q1 de 
dimensions finies, il faut rajouter les conditions de transversalité : 


p(0) —+— Tao) Qo 
PT) L Tir) Qi 


L'idée sous-jacente est que, si le point final ou de départ est libre, alors parmi tous les pro- 
blèmes point-point possibles, il faut choisir celui (ou ceux) qui minimise(nt) le coût. 


À titre d'exemple, nous allons appliquer le PMP sur un cas simple. Nous étudierons les tra- 
jectoires du système de Dubins (2.2) qui minimisent la fonctionnelle J représentant un compromis 
pondéré entre la longueur et l’énergie du signal de contrôle u(-) d’une trajectoire q : [0, T] — Q. 
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Celle-ci est définie par 
T 
Hu) | 1 + au(t)?dt, 
0 


avec & > 0 une constante de pondération. 
Les trajectoires minimisant un critère similaire ont été étudiées par Yuri Sachkov et al. dans une 
série d'articles [92, 116, 117] pour un système de type Reeds-Shepp. 


Écrivons le problème (Pz,,), pour le système de Dubins, associé au coût L,(u) = 1 + au? : 


(Pz,,) Minimisation d’un coût intégral 


par rapport à tous les contrôles admissibles u(-) € L°([0,T],7) conduisant le système (2.2) 
d’un point de départ go à un point d’arrivée gap en un temps T libre. 


Soit le Hamiltonien 
HN P, qu) = px COS(0) + pysin(0) + pau + AL (u), (2.9) 
avec p = (Pr, Py, Do) € R° et À < 0. 


Si (À,p(-),q(-),u(-)) une extrémale du problème (P7,,) alors elle vérifie les conditions du PMP, 
dont le système (2.7). 


Comme le temps final Test libre, le Hamiltonien est identiquement nul le long de la trajectoire 
optimale [3], i.e. 
0 = p4 cos 0 + p, sin 0 + pou + ALA(u). (2.10) 
L'équation de q dans (2.7) n’est rien d’autre que (2.2), quand à l’équation de p, aussi appelée 
équation adjointe, elle peut explicitement être écrite 


Py = 0 11) 
Do = Px Sin 0 — p,, cos Ô 


Remarque 2.15. Par la suite, les valeurs constantes de p,(-) et p,(-) seront notées, respectivement, 
Pr et Py- 


En utilisant les équations (2.2), (2.8), (2.10) et (2.11) nous sélectionnons les trajectoires candi- 
dates à être solution de (Pr, ). 


Dans le cas où une synthèse complète du problème considéré a été obtenue (comme par exemple 
dans l’article [117] pour le système de Reeds-Shepp avec un critère similaire à celui présenté précé- 
demment), il suffit d'appliquer les règles obtenues lors de la synthèse pour construire la trajectoire 
optimale. 
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FIGURE 2.3 — Résultats de planification optimale, utilisant le PMP, minimisant le compromis L, en 
temps libre (la trajectoire continue correspond à a& = 10 et la trajectoire en pointillé à & = 2). Le 
point initial est go — (9,2,0) et le point final est qcible — (5,6, 7/2). 


Si ce n’est pas le cas, il faut calculer la synthèse numériquement en intégrant le système hamil- 
tonien (2.7). La trajectoire solution est construite en cherchant la valeur du vecteur adjoint initial 
(dans notre cas p(0) = (p+, py, po(0))) et du temps final T, s’il est libre. Si le temps est fixé, l'équation 
(2.10) n’est plus valable. 

Les algorithmes que nous utilisons pour obtenir numériquement la synthèse sont des algorithmes 
de tir. Pour plus d'informations, voir [121, 135, 93, 83]. 


Deux solutions du problème (P7,,) sont représentés en Figure 2.3. 

Dans ces expériences, deux points arbitrairement fixés sont reliés par une trajectoire. Le temps 
est libre et le contrôle du cap n’est pas contraint, i.e. u € R. La différence entre les expériences est la 
valeur de pondération du compromis utilisée. Nous constatons que pour obtenir la trajectoire conti- 
nue, où la valeur de a est la plus grande, l’énergie du contrôle est moins importante (la trajectoire 
à une courbure moins importante). 


Beaucoup de travaux ont été effectués pour produire des synthèses de problèmes optimaux mini- 
misants le temps. En premier lieu, nous pouvons citer le travail de U. Boscain et B. Piccoli qui, dans 
leur livre [31], exposent une méthode de construction des synthèses optimales en temps pour des 
systèmes de dimension 2. Cet ouvrage a comme point de départ les résultats de [32, 33, 103, 125]. 

Les premiers travaux sur les trajectoires temps-minimal pour les systèmes de Dubins et de Reeds- 
Shepp ont été effectués, respectivement, par L.E. Dubins [51] et Reeds et Shepp [110]. Ils ont carac- 
térisé les trajectoires optimales en temps avec un contrôle borné. 

Les trajectoires temps-minimal du système de Reeds-Shepp ont été étudiées par Sussmann et 
Tang [127] ainsi que Boissonnat et al. [27]. Souères et Laumond ont construit la synthèse temps- 
minimal de ce système [122]. 
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La synthèse temps-minimal pour le système de Dubins a été obtenue par Bui, Boissonnat, Souères 
et Laumond [31|. Souères et al. étudièrent aussi les trajectoires optimales en temps pour rejoindre 
une droite et ont produit la synthèse associée [14]. Dans le même temps, Bakolas donna une synthèse 
temps-minimal pour le système de Dubins en présence de vent [13]. 

Bien d’autres problèmes de planification ont été étudiés grâce à ce principe [35]. 


Dans le paragraphe 3.2, nous allons calculer la synthèse du problème temps-minimal pour re- 
joindre un cercle, pour le système de Dubins. 


2.3 Le suivi de convoi 


Dans cette partie, nous nous intéressons au problème de planification de trajectoires dans le cas 
d’un suivi de cible mobile. La problématique est donc d'utiliser les méthodes présentées précédem- 
ment pour prendre en compte la non stationnarité du point final. 

La planification de trajectoires vers un point final en mouvement permet de répondre à certains 
besoins spécifiques des missions effectuées par un drone (e.g., espionnage, protection). Pour cela nous 
pouvons aussi utiliser les patterns définis par la norme mise en place par l’organisation de l’aviation 
civile internationale [95]. 


Nous ne nous intéressons qu’au contrôleur de haut niveau du vecteur. En d’autres termes, nous 
supposons qu'il existe un module embarqué qui permet au vecteur de suivre une trajectoire prédéfi- 
nie (par points de passage par exemple). C’est le cas dans [108, 25] où la trajectoire fournie par un 
module de planification est transmise à un contrôleur de vol du commerce, qui effectue le contrôle 
de bas niveau (cf. Figure 1.1). 


Différentes méthodes de suivi de cibles sont décrites dans [143, 42]. 
Classiquement, les auteurs utilisent des méthodes stochastiques [22], des méthodes d'optimisation 
numérique [78] ou encore des méthodes basées sur un modèle de type proie/prédateur [12]. 


2.3.1 Méthode quadratique 


Dans [29], les auteurs présentent une méthode quadratique de contrôle qui permet de prendre 
en compte les contraintes sur l’état et la commande. Nous l’adaptons à la planification de trajectoires. 


Nous supposons que le système est affine dans le contrôle, 1.e. il est de la forme (2.3). La méthode 
consiste à trouver un contrôle minimisant l'écart entre la vitesse q et la valeur cible cible. 

Il faut alors chercher à exprimer cet écart en fonction de la commande u. Le but est de résoudre 
le problème de minimisation quadratique associé, avec ou sans contraintes sur l’état et/ou sur la 
commande [105]. 


La commande u*(-) solution du problème de minimisation permettant de construire la trajectoire 
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peut répondre, par exemple, au problème de la forme (cf. paragraphe 3.1) 


(Z- Re) A(Z Ru‘) = min L(Z — Ru) A(Z — Ru)|. 


où Z, R sont des matrices et À est une matrice de pondération. 


Ce problème de minimisation est standard et peut être résolu très simplement grâce à la pro- 
grammation quadratique (e.g., par une méthode de relaxation [105]). 


Dans le paragraphe 3.1, nous utiliserons cette méthode quadratique pour répondre au problème 
de couplage vecteur/capteur. 


2.3.2 L'utilisation de primitives 


Le problème de planification, en temps réel, est soumis à de nombreuses contraintes, dont celle 
du temps de calcul. C’est pourquoi, beaucoup de modules de planification sont basés sur des trajec- 
toires types, ou primitives, que peut effectuer le vecteur. La trajectoire planifiée est alors composée 
d’une concaténation de trajectoires types. 


Ces primitives peuvent être construites à partir de la résolution d’un problème de contrôle op- 
timal. Par exemple, Latombe a utilisé le modèle de Reeds-Shepp, et ses courbes temps-minimal, 
pour construire un algorithme de planification rapide pour une voiture autonome, dans un milieu 
encombré [81]. 

Frazzoli et al. ont proposé des trajectoires types, appelées trim primitives, pour un hélicoptère 
à 6 degrés de libertés. À partir de ces primitives de mouvement, obtenues par le biais d’un module 
de pilotage automatique, ils construisent un module de planification concaténant ces morceaux de 
trajectoires pour optimiser un certain coût [58]. 

Cette approche a été utilisée pour la planification de trajectoires de robots bipèdes. Par exemple, 
Hauser et al. proposent une méthode probabiliste de planification utilisant des primitives de mou- 
vement d’un robot humanoïde [68]. 

Enfin, dans [85, 70], les auteurs utilisent des primitives de trajectoires définies à partir de la 
différence de vitesse entre le vecteur et sa cible. 


2.3.3 Des méthodes basées sur des considérations géométriques 


Pour certaines missions, des patterns prédéfinis sont proposés et, dans ce cas, des algorithmes de 
planification temps-réel peuvent être mis en œuvre. 


Une première méthode, décrite dans Rañi et al. [107], consiste à donner un cap à suivre au vecteur 
permettant de rejoindre un cercle (à la même altitude que le vecteur) centrée sur la cible considérée. 
C’est la méthode de la tangente qui n’est utilisable que si le vecteur est à l’extérieur du cercle visé. 
Dans le cas où le vecteur est à l’intérieur du cercle, une autre méthode de planification est appliquée, 
par exemple une méthode de stabilisation basée sur le principe de LaSalle [37]. 
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FIGURE 2.4 - Schéma du principe de la méthode de la tangente 


Lorsque le vecteur est en dehors du cercle, l’idée est de considérer les deux tangentes au cercle, 
passant par le centre de masse du drone. Le cap Ouav.. que doit suivre le vecteur est celui qui est le 
plus proche du cap actuel et qui permet de suivre une des deux tangentes. En utilisant les notations 
de la Figure 2.4, Ouavanre = Ouav + Min(d, 042), avec 0x1, 04 les angles entre le cap actuel et les 
tangentes au cercle. Cette commande de cap est envoyée au module de pilotage automatique du 
drone. 


Une autre méthode, basée sur des considérations géométriques, est développée dans l’article 
de Theodorakopoulos et al. [130]. La problématique concerne les drones HALE et le couplage vec- 
teur/capteur. Le vecteur doit rejoindre la projection du point cible dans son plan d’altitude constante 
avec une trajectoire maximisant le temps d’observation de la cible. La caméra embarquée est soumise 
à des contraintes sur son champ de vision. 


Le cap que doit prendre le drone est calculé de sorte que : 


e La cible reste au centre du champ de vision de la caméra, 
e La distance entre le vecteur et la cible soit plus grande que le rayon de courbure minimal. 


La loi de guidage latéral permet de contrôler le cap du vecteur. 


Après une étude géométrique du problème, les auteurs constatent qu’une trajectoire en spirale 
permet d'obtenir un temps de visibilité de la cible maximal. 
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2.3.4 Les méthodes utilisant le principe de LaSalle 


Comme mentionné au paragraphe précédent, certains auteurs utilisent des méthodes basées sur 
le principe de LaSalle pour stabiliser le drone sur un pattern. 


Wood, dans sa thèse [142], ainsi que Hamel et al. , dans [63, 89, 71, 21], ont considéré une méthode 
de suivi de trajectoire, basée sur le principe de LaSalle, pour un véhicule à décollage vertical. 

Lawrence et al. ont appliqué de telles techniques sur des modèles de drone de type HALE 
[84, 59, 60]. 


Ce principe de LaSalle est aussi utilisé pour trouver des méthodes de suivi de cible décrivant des 
trajectoires types (e.g., des droites ou des cercles) [100, 144]. 


Dans tous les travaux cités, les méthodes de stabilisation sont appliquées en supposant que le 
problème de suivi concerne la stabilisation dans le plan d’altitude constante du drone. Nous allons, 
dans le paragraphe 3.2, proposer une méthode, basée sur le principe de LaSalle, permettant de passer 
outre cette supposition et susceptible d’être appliquée à des drones de type MALE. 


2.3.5 La résolution du problème de suivi via le PMP 


Le Principe du Maximum de Pontryagin a déjà été exploité pour la mise en œuvre d’algorithmes 
de suivi de convoi. 

Par exemple, les trajectoires en temps-minimal, pour le système de Dubins, sont utilisées par 
Bhatia et al. [25, 24]. Dans ces articles, les auteurs considèrent plusieurs drones et veulent qu’ils 
rejoignent le même point, au même instant. Cet exemple peut être étendu, par la suite, pour effectuer 
de la reconnaissance en coopération. 

Ce même problème a été étudié, plus récemment, par Ortiz et al. dans [98], pour faire de la 
surveillance de zone à plusieurs drones. 


Ding et al. ont considéré le problème de la protection d’un convoi en utilisant des trajectoires 
temps-minimal, du système de Dubins, rejoignant un cercle centré sur le convoi [8]. Dans cet article, 
le rayon de courbure minimal des trajectoires du drone est supposé plus grand que le rayon du cercle 
de protection autour de la cible si bien que le vecteur ne peut pas suivre le pattern autour du convoi. 


Dans le paragraphe 3.2, nous proposerons une synthèse des trajectoires temps-minimal, pour un 
drone de type HALE, permettant de rejoindre le cercle de rayon le rayon de courbure minimal du 
vecteur. Cette synthèse peut être utilisée pour la protection des convois. 


2.4 Prédiction du mouvement de la cible 


Dans la partie précédente nous avons proposé quelques méthodes pour rejoindre une cible. Dans 
le cas d’une mission de protection d’un convoi (une cible amie), la position cible qi est considérée 
connue (le convoi peut communiquer ses coordonnées au drone). Cependant, dans le cas d’une cible 
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ennemie, les méthodes ont besoin d’un module supplémentaire qui permet d’obtenir les informations 
nécessaires. 

La plupart des auteurs, pour estimer la position de la cible, utilisent des méthodes probabilistes 
[42, 86]. 


Dans cette partie, nous présentons, succinctement, deux méthodes couramment utilisées. 


2.4.1 Le filtre de Kalman 


Une méthode classique d’estimation de la position future de la cible est basée sur un filtre de 
Kalman étendu. Ce filtre estime les états d’un système dynamique à partir d’une série de mesures 
incomplètes et bruitées. Le filtre de Kalman étendu est utilisé pour observer des systèmes non- 
linéaires : c’est une extension du filtre de Kalman linéaire basée sur la linéarisation autour de la 
trajectoire estimée. 


Ce filtre est présenté, par exemple, dans le livre de Besançon [23] et est appliqué sur un quadcopter 
dans l’article de Boizot et al. [120]. 


Remarque 2.16. Le filtre de Kalman étendu ne converge pas globalement. Il existe, par contre, des 
preuves de convergence locale [15]. 


Cette méthode de prédiction de la position d’une cible est utilisée, par exemple, dans l’article de 
Prévost et al. [106]. Nous l’utilisons dans le paragraphe 3.2 pour réaliser du suivi de convoi. 


2.4.2 Le filtrage particulaire 


Le filtrage particulaire, ou méthode de Monte-Carlo séquentielle, est une technique d’estimation 
de modèles basée sur des simulations probabilistes. 


Cette méthode est détaillée dans [83, 132, 86, 129]. Les principales étapes sont les suivantes : 


1. Initialisation : fixer un certains nombres de particules p, correspondant à des états initiaux 
distribués uniformément dans l’espace d’état et calculer leurs poids associés (en fonction de 
leurs vraisemblances), 


2. Propagation : calculer les mouvements de chaque particule en fonction d’un modèle probabiliste 
[19], puis calculer leurs poids associés, 


3. Sélection : tirer avec remise p particules, dans l’ensemble des particules existantes, proportion- 
nellement à leurs poids. Associer le poids 1/p aux p nouvelles particules. Retourner à l’étape 
2. 


Le modèle de mouvement Brownien est un exemple de modèle appliqué à l’étape 2. 


Le filtrage particulaire a été utilisé par Hongda et al. dans [37] pour estimer la position d’une 
cible mobile. 
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2.5 L’évitement d’obstacles 


Toutes les méthodes proposées jusqu’à présent ont été développées dans le but d’être utilisées 
dans un environnement non obstrué, i.e. dans le cas où le vecteur peut se déplacer dans tout l’espace. 

Dans la réalité, l’environnement est souvent obstrué. Dans ce paragraphe, nous verrons comment 
détecter des obstacles puis planifier la trajectoire dans un environnement contraint. 

Le terme obstacle désigne, ici, aussi bien des entités physiques (e.g., bâtiments, reliefs, …) que 
virtuelles (e.g., zones de non vol, frontière, ….). 


Afin d'éviter ces obstacles, il est nécessaire de les localiser. Si l'emplacement des obstacles virtuels 
est supposé connu, ce n’est pas le cas pour les obstacles physiques. En effet, à tout moment, un objet 
peut entrer dans l’environnement proche du vecteur. 

Il est important de modéliser ces obstacles. Nous souhaitons créer une entité permettant de 
définir un obstacle (physique ou virtuel) pouvant être reconnu par le module de planification. 

Enfin, une fois les obstacles définis, il faut utiliser les méthodes de planification connues afin de 
proposer une trajectoire réalisant nos objectifs (e.g., les objectifs du couplage vecteur/capteur) et 
évitant d'entrer en collision avec les obstacles. 

Ces trois points seront traités, succinctement, dans cette partie. Pour une bibliographie plus 
complète, se référer à [83]. 


2.5.1 La détection des obstacles 


Nous allons ici énoncer quelques méthodes, couramment utilisées, pour construire un module de 
détection d'obstacles physiques. 


Le plus naturellement possible, pour mimer le comportement humain, la vision peut être exploitée 
pour construire un module de détection. En effet, dans la vie courante, l’homme utilise principalement 
sa vision pour se repérer et repérer les obstacles proches. 

Kim et al. font de la reconnaissance de scène grâce à de la vision. Dans une série d’articles 
[79, 109, 59], ils ont mis en place, sur un drone, un module de détection et de suivi de routes via 
une caméra embarquée. 

Par ailleurs, Kim et al. ont implanté sur un drone un module de recherche et de suivi de rivière, 
basé sur du traitement d'images [108]. 


Pour obtenir de meilleurs résultats, la caméra peut être couplée à un radar. Ce couplage de 
capteurs est appliqué à la détection d’objets et à la reconstruction de scènes [124]. 


Au lieu d’un radar, certains travaux utilisent un laser pour permettre une détection sur un envi- 
ronnement plus large. C’est le cas, par exemple, dans [101] où l’auteur présente une caméra couplée 
à un dispositif laser pour effectuer de la détection d'obstacles, pour un problème d'assistance à la 
conduite. 
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FIGURE 2.5 — Schéma d’une méthode simple d’évitement d'obstacles. La trajectoire relie les points 
go et q1 et évite les obstacles (en noir). 


Nous ne développerons pas plus ce paragraphe puisque, par la suite, nous supposerons que la 
position des obstacles est connue. Pour plus de détails sur le sujet, se référer à [64]. 


2.5.2 La discrétisation de l’espace 


En considérant la position des obstacles connue (physiques ou virtuels), il faut les modéliser et 
les prendre en compte dans l’algorithme de planification. 


En appliquant une méthode simple de planification, comme la méthode de la tangente, les obs- 
tacles sont modélisables par des cercles [37]. Ces cercles doivent englober les obstacles et avoir un 
rayon plus grand que le rayon de courbure minimal du drone. Le cap objectif à envoyer au module 
de pilotage automatique est celui qui permet de suivre les tangentes aux obstacles situés entre le 
point cible et le vecteur (courbe continue sur la Figure 2.5). 


Le problème principal de cette méthode est le suivant : si l’environnement du drone est constitué 
d’un couloir étroit, alors nous aboutissons à l’obstruction du couloir. Par exemple, sur la Figure 2.6 
l'algorithme de planification considérera que le drone ne peut pas passer entre les deux obstacles 
puisque leurs cercles circonscrits se chevauchent au niveau du couloir. 

Ce problème d’obstruction peut aussi survenir en considérant d’autres formes géométriques 
simples pour modéliser les obstacles (e.g., des carrés ou des triangles). 


Une méthode plus efficace de modélisation des obstacles passe par une discrétisation plus fine 
de l’espace. L'idée est alors de considérer des formes géométriques, que nous appelons cellules, per- 
mettant de découper l’environnement. Les obstacles sont alors constitués d’un ensemble de cellules. 
Une discrétisation simple peut être effectuée par des carrés, e.g., les pixels de l’image de la Figure 2.7. 
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FIGURE 2.6 — Schéma d’une configuration d'obstacles ne convenant pas à la modélisation par des 
formes géométriques simples. 


La décomposition par cellules identiques, i.e. l’environnement est discrétisé uniformément par 
une forme simple, peut poser problème. En effet, il peut y avoir une majorité d’éléments ne conte- 
nant pas d'informations utiles (e.g., les pixels blancs de la Figure 2.7). Une idée serait alors de les 
regrouper pour constituer des cellules plus grandes (pas nécessairement identiques) pour modéliser 
l’espace. La décomposition de Boustrophédon [50, 91] le permet. 


Cependant, il existe une méthode de discrétisation plus efficace. Celle-ci est basée sur la trian- 
gulation de Delaunay et les graphes de Voronoï [137, 96]. Elle utilise la théorie de Morse discrète 
introduite dans [57]. L'idée est d'utiliser les graphes de Voronoï afin de créer une cellule par obstacle, 
la frontière de ces cellules étant toujours équidistante de tous les obstacles. Une stratégie de guidage 
sécuritaire est de suivre cette frontière. Un exemple d'utilisation de cette discrétisation est présenté 
en Figure 2.8. 

Cette triangulation a été utilisée, par exemple, par Belta et al. dans [18]. Nous l’utiliserons dans 
la prochaine partie ainsi que dans le paragraphe 3.4. 


2.5.3 Les méthodes de planification en milieu contraint 


Une fois une modélisation des obstacles obtenue, par discrétisation de l’espace par exemple, nous 
appliquons une méthode de planification. 


Généralement, après cette discrétisation de l’espace, des algorithmes de planification discrets 
sont appliqués. La majeure partie de ceux-ci sont basés sur la théorie des graphes [73]. 

Certains auteurs planifient le mouvement du vecteur en combinant les éléments d’un alphabet 
de trajectoires [18]. D’autres utilisent des algorithmes de cheminement, plus communément appelés 
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FIGURE 2.7 — Schéma d’un environnement obstrué, discrétisé via une pixelisation. 


FIGURE 2.8 — Exemple d’un graphe de Voronoï. 
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FIGURE 2.9 —- Exemple d’un fil d’Arianne (en pointillés bleu). 


algorithmes de déménageur de piano [50, 73], tels que l’algorithme de Dikjstra [1, 136] ou le A* 
[66, 140]. 

Un exemple d'utilisation de l’algorithme A* est proposé dans [69]. Des points de passage sont 
calculés via cet algorithme et sont fournis au contrôleur bas niveau du drone, basé dans ce cas sur 
un PID. Un exemple de trajectoire fournie par cet algorithme est représenté sur la Figure 2.9. Sur 
cette figure, le résultat de la planification est représenté en pointillés. La discrétisation est construite 
à partir du graphe de Voronoï, cf. Figure 2.8. 


La méthode présentée précédemment, utilisant des algorithmes de cheminement, permet de cal- 
culer une trajectoire évitant les obstacles. Cependant les trajectoires calculées ne sont pas toujours 
admissibles pour le vecteur car sa dynamique ne permet pas de les suivre. 

Une idée serait alors d'appliquer le PMP avec contraintes sur l’état, sur le système cinématique 
du drone [135, 43, 104]. Cette méthode permettrait d'obtenir une trajectoire pouvant être suivie par 
le drone. De plus, en considérant des contraintes sur l’état, i.e. des ensembles auxquels l’état ne doit 
pas appartenir, la trajectoire évitera aussi les obstacles. 


Bien que l’idée d'obtenir une trajectoire optimale soit intéressante, le PMP contraint est assez 
difficile à mettre en œuvre. Afin d'utiliser les avantages du PMP (l’optimalité des trajectoires au 
sens d’un certain coût et l’utilisation de contraintes dynamiques), une méthode a été mise en œuvre 
par Latombe dans [81]. La planification se fait en deux étapes : 


1. La première étape consiste à calculer une trajectoire, grâce à un algorithme de cheminement, 
reliant les points de départ et d’arrivée, tout en évitant les obstacles. Nous nommerons ce 
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FIGURE 2.10 — Principe de l’algorithme d’évitement d'obstacles 


chemin le fil d'Ariane qui peut être, par exemple, le résultat de l’algorithme A* en Figure 2.9. 


2. La deuxième étape consiste à utiliser ce fil d'Ariane comme support pour calculer une tra- 
jectoire admissible pour le véhicule en utilisant, par exemple, le PMP ou des primitives de 
trajectoires optimales en temps. 


En résumé, l’idée sous-jacente est de réduire le problème contraint à des sous-problèmes sans 
contrainte sur l’état. 


Nous avons utilisé cette idée pour construire un algorithme d’évitement d’obstacles. 
Tout d’abord, cet algorithme doit trouver des points de passages dans le but d’appliquer une 
méthode de planification sans contrainte et ainsi obtenir une trajectoire admissible. La recherche de 
ces points se fait grâce à l’algorithme discret A* qui permet de trouver un fil d'Ariane, reliant les 
points de départ et d’arrivée, et respectant les contraintes d'état. 
Le PMP est alors appliqué pour trouver une trajectoire optimale, en fonction d’un coût choisi 
par l'utilisateur, reliant la position courante du vecteur à un point du fil d'Ariane. 
Voici les étapes de construction de la trajectoire. Celle-ci est illustrée sur la Figure 2.10. 
Étape 1 Initialisation du point courant (la position du vecteur), 
Étape 2 Recherche du fil d'Ariane, par discrétisation de l’espace et en utilisant l'algorithme A* 
(courbes annotées 2 et 2’ sur la Figure 2.10), 

Étape 3 Recherche d’une trajectoire, avec le PMP, reliant le point courant et un point du 
fil d'Ariane, le plus éloigné possible. Les courbes annotées 3, 3’ et 4, 4’ de la Figure 2.10 
représentent cette étape e 

Étape 4 Mise à jour du point courant et retour à l'étape 2, tant que celui-ci n’est pas proche 

du point d’arrivée (sur la Figure 2.10, cette boucle est effectuée deux fois). 


1. Les courbes annotées 3, 3° correspondent à un prétraitement effectué localement pour améliorer la recherche des 
courbes optimales (annotée 4, 4’). 
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FIGURE 2.11 — Recherche d’une trajectoire optimale, en terme de longueur et de courbure, avec un 
temps final libre dans un espace contraint. 


Cette méthode de recherche de trajectoires à été implantée sous Matlab. La Figure 2.11 corres- 
pond à une planification vers un point final fixe, en temps libre. La courbe en pointillés représente 
la trajectoire non-admissible utilisée comme fil d'Ariane pour calculer une trajectoire admissible et 
évitant les obstacles (courbe continue). 
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CHAPITRE 3. RÉSUMÉS DES CONTRIBUTIONS PERSONNELLES 


3.1. Comment effectuer le couplage vecteur/capteur ? 


Jusqu'ici nous n’avons présenté qu’un panorama des méthodes de planification de trajectoires. 
Nous traitons maintenant le problème du couplage vecteur/capteur. 

Au paragraphe 3.1, nous montrons que le problème du couplage vecteur/capteur peut être résolu 
simplement en cas de contraintes en position sur la caméra. Dans cette même partie, nous proposons 
une méthode permettant de prendre en compte la caméra lorsque celle-ci n’est soumise à aucune 
contrainte de mouvement ce qui nous amènera à constater que le problème du couplage se situe 
essentiellement dans la planification du vecteur (pour éviter de perdre de vue la cible à cause d’un 
mauvais positionnement). 

Les questions qui se dégagent sont les suivantes : 


1. Comment effectuer le couplage vecteur/capteur ? 


2. Comment définir une trajectoire permettant de rejoindre un pattern nous positionnant dans 
des conditions optimales pour l’observation de la cible ? 


3. Comment suivre une cible mobile pour l’observer au mieux ? 
4. Comment suivre une cible mobile en milieu contraint (e.g., en présence d’obstacles) ? 


Dans la suite, nous donnons une réponse à ces problématiques. Les réponses aux questions 2 et 3 
sont détaillées dans les articles présentés au paragraphe 3.2 et dans les Annexes À et B. Ces articles 
sont publiés ou en cours de publication. 


3.1 Comment effectuer le couplage vecteur/capteur ? 


Dans ce paragraphe, nous proposons une stratégie de couplage vecteur/capteur dans le cas d’un 
point final fixe ou mobile. 
Cette stratégie de contrôle est présentée au paragraphe 3.1.1. Elle induit une première pré-étude 
du problème de couplage vecteur/capteur (cf. paragraphe 3.1.2) : 
e Prise en compte des contraintes en position de la caméra pour effectuer la planification du 
vecteur. 
e Dans le cas où la caméra n’est pas contrainte, le problème de couplage revient à étudier 
uniquement la planification du vecteur (cf. paragraphe 3.2). Dans ce cas, nous suggérons une 
méthode indépendante de contrôle de la caméra (cf. paragraphe 3.1.3). 


3.1.1 Description de la méthode 


L’algorithme utilisé pour répondre au problème de planification posé au paragraphe 1.1 est 
inspirée de [29] (introduit au paragraphe 2.3.1). Celui-ci vise à trouver une trajectoire qui permet 
de prendre en compte les contraintes sur le vecteur et la charge utile. 

Par la suite, nous considérons un drone HALE qui a une caméra pour charge utile. Le cas des 
drones MALE est similaire. 

Nous supposons connu l’état de la cible Giple = (cible: Ycible), ainsi que ses dérivées dcible €t Qcible: 
Ces informations sont données par la cible elle-même (dans le cas du suivi de convoi ami), ou bien 
par un module d’estimation, cf. paragraphe 2.4. 
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3.1. Comment effectuer le couplage vecteur/capteur ? 


L'idée est de considérer le problème comme une minimisation de l’erreur entre la vitesse courante 
de l'UAV et la vitesse actuellement connue de la cible cible (e.g., un convoi à suivre). 


Nous utilisons un modèle de Dubins légèrement modifié pour décrire la cinématique du drone 
HALE et la dynamique de sa caméra (dans le plan). 


& = vcos 

ÿ = vsin 0 

ê=u (3.1) 
v=p 

à=u+ ue +(a Ô) 


avec q = (x,y,0,v,a) l’état du système et U = (u,u,u.) € # ses commandes. Nous posons # — 
este) X brie fase] X le, Gras] AVEC max; Umax POSitiFS. 

Remarquons ici que la vitesse du drone v peut évoluer grâce à la commande d'accélération y. 
Nous la supposons contrainte et vérifiant 0 < Unin < U < Umax- 

Comme pour le système de Dubins classique, le point (x,y) € IR? définit les coordonnées du drone 
dans le plan d'altitude constant et 0 son cap. Les commandes 1 et uw sont associées au vecteur et 
régissent, respectivement, son accélération et sa commande de lacet. Un schéma est donné en Figure 
3.1. 

L’angle à régit le positionnement de la caméra dans le plan. Nous considérons que cet angle 
définit l’angle de la caméra dans le repère global, en d’autres termes a = & + 4 où & est l’angle de 
la caméra dans le repère associé au drone. La commande u. régit la position de la caméra dans ce 
repère, i.e. u. commande &. Nous supposons que la caméra est un système linéaire du premier ordre 
(la caméra est dirigée par un moteur électrique qui peut être considéré comme tel) : 


a 4h L 
À = —Ue — —À, 
Te Te 
avec K le gain statique du système et 7. sa constante de temps. Dans la suite, nous supposons que 
K = 1 ce qui peut toujours être effectué quitte à changer les bornes du contrôle uw. 
Remarque 3.1. L’équation de la caméra étant une équation différentielle du premier ordre, les 
contraintes de position sur & sont celles sur la commande u,, i.e. nous supposons &min < 4 < Amax- 


Grâce aux connaissances sur la cible, nous calculons œible €t cible, OÙ &cibie est l’angle que doit 
avoir la caméra pour observer la cible. 


: Ycible — Y 
cible = AICtan | 
Tcible — L 


Fr (Habite Po y )(Sébie _ #) _. (Ysibie = y)(£cible _ ä) 

ces (Laibie — 6 ad (cite — y}? 

Remarque 3.2. Les cas particuliers où æcjple — 2 = 0 où (teible — 7)? + (Ycible — Y)? = 0 sont traités 
à part. Dans le premier cas nous posons ace — T. Lorsque (æcible — 7)? + (Ycible — y)? = 0 nous 
considérons que &cible = 0. 
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3.1. Comment effectuer le couplage vecteur/capteur ? 


FIGURE 3.1 — Schéma du drone HALE et de la charge utile considérée 


Nous minimisons les erreurs de positionnement € et €. du vecteur et de la caméra, par rapport 
aux objectifs voulus. Ces erreurs sont définies de la manière suivante : 


.: Tcible — À 
Ucible — Y 
Ece — cible — 


Le système régissant la position du vecteur (3.1) est invariant sous l’action du groupe des dépla- 
cements du plan, SE(2) [8]. Nous supposons donc, en utilisant la méthode décrite par P. Rouchon 
[113], que l'erreur dynamique de suivi de cible peut s'exprimer comme un système linéaire d'ordre 
deux, i.e. l'erreur € vérifie 

Ë = LÉ + 12€, 


avec L1 et 2 deux constantes réelles. 

En ce qui concerne l’erreur de suivi de la cible par la caméra, sachant que celle-ci est repré- 
sentée par un système linéaire d’ordre 1, nous considérons que l'erreur est solution d’une équation 
différentielle linéaire d'ordre 1, i.e. l’erreur €, vérifie 


Ée — L3€c; 


avec L3 une constante réelle. 
Les deux conditions sur l’évolution des erreurs conduisent au système suivant 


O0 = Ê — spé — 19€ 
0 = ë, — 13€ 
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3.1. Comment effectuer le couplage vecteur/capteur ? 
La solution stationnaire (£,£c) = (0,0) est stable si s3 < 0 et les valeurs propres de la matrice 
0 1 
Lo 11 


sont à parties réelles négatives. 
De plus, en utilisant (3.1), nous avons 


Ec — cible — 


. ; 1 
Ec — cible — U — 0. + —(a = 0) 


C C 
_ Li _ () 
Ucible y 
: : ( COS à 
Ycible v sin Ô 
É cible —v sin Ü cos Ô 
— | . u— |. nl 
Ucible v cos Ô sin 0 


En résumé, les commandes u, 1 et u. sont solutions du système de contrôle 


et 


Tcible Tcible — v cos Ô Tcible — L ( 
0 — Ücible — 11 | Ycible — VSINO | — 19 | Ycible — Y | — 43 0 
: 1 
cible + (a . 0) (ù 0 cible — 
—v sin Ô cos 0 0 
— vcos0 | u— | sing | — | O | u,, 
1 
1 0 _ 
qui se réécrit sous la forme 
0 = Z — RU, (3.2) 
où 
Tcible Tcible — V COS Ü Tcible — TL 0 
Z = ee Uobe—v8m0 | — | - 9 | 8 0 ; 
cible + 0 0 cible — 
—vsin 6 cosû 0 
R= | vcosô sinÿ 01, 
1 0 + 
ns 
U = (u, H; üc). 


La commande U*(-) solution du système (3.2) vérifie aussi le problème de minimisation quadra- 


tique suivant, pour tout & > 0 


(Z() - RHU*()" A (Z(6) — ROU*(E)) = qu (ZE) — RU)" A(Z(6) - ROU)|, 
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où À est une matrice de pondération que nous utilisons pour pondérer l'erreur de positionnement 
du vecteur par rapport à celle de la caméra, ou inversement. 


Ce problème de minimisation a une solution explicite triviale si U n’est pas contraint. Dans le 
cas contraire, nous sommes dans le cadre de la programmation quadratique (e.g., la méthode simple 
de relaxation peut être appliquée [105]). 

Afin d’asservir la vitesse v, que nous ne commandons pas directement, nous modifions les bornes 
de la commande d’accélération 4 en fonction de la vitesse en cours. Ces deux bornes définissent donc 
deux fonctions à un paramètre telles que 


HUE — max tanh(o(v + Umin)) 


U+ UE. LUimax taNh(o(v — Umax)), 


avec © une constante réelle. 


3.1.2 Application de la stratégie de couplage vecteur/capteur 


La méthode quadratique de recherche de trajectoires a été implantée sous Matlab. Les Figures 
3.2, 3.3, 3.4, 3.5 sont des exemples de résultats de suivi de cibles. Elles correspondent à des contraintes 
différentes sur la caméra, ainsi qu’à des paramètres de couplage vecteur/capteur différents. C’est à 
dire que la matrice de pondération À est différente et le champ de vision de la caméra est contraint 
ou non. 


Sur chacune des Figures 3.2 à 3.6, la sous figure de gauche représente la trajectoire du vecteur 
obtenue par la méthode quadratique et la sous figure de droite représente l’erreur de positionnement 
de la caméra en fonction du temps. Sur chacune des sous figures de gauche, la trajectoire en trait 
continue est celle du vecteur et la trajectoire en traits pointillés correspond au mouvement de la 
cible à suivre. Les moyennes et écart-type des erreurs de positionnement de la caméra (& et 04, 
respectivement) et du vecteur (|| et le; respectivement) sont indiqués pour chaque simulation. 


Les Figures 3.2 et 3.4 correspondent à des simulations avec une matrice de pondération qui donne 
une importance plus grande au positionnement de la caméra. À l'inverse, pour obtenir les Figures 
3.3 et 3.5, la position du vecteur est jugée plus importante. 

Les deux premières figures (les Figures 3.2 et 3.3) sont des simulations de planifications avec des 
contraintes sur le positionnement de la caméra. Ces contraintes n’ont pas été appliquées pendant les 
simulations correspondant aux Figures 3.4 et 3.5. 


Un dernier exemple d’application est présenté sur la Figure 3.6. Lors de cette simulation, nous 
avons utilisé cette méthode pour la mission suivante : 
e Le vecteur ne doit pas franchir la frontière (matérialisée par une ligne en tiret-point). 
e La caméra doit suivre le pattern représenté par des tirets. L’angle d'observation de la caméra 
est supposé contraint : à € [—-1/3,7/3]. 


46 


CHAPITRE 3. RÉSUMÉS DES CONTRIBUTIONS PERSONNELLES 


3.1. Comment effectuer le couplage vecteur/capteur ? 


D 


So a 
T 
Rs 
= 
nes 


_2t 


0 5 10 15 20 25 0 5 10 15 20 25 30 35 20 25 50 
T Temps 
(a) [ll = 2.021, oyey = 12077 (b) & = 0.2226, 0. = 1.1529 


FIGURE 3.2 — Trajectoire primant l’erreur de positionnement de la caméra. Le champ de vision de 
la caméra est supposé contraint : à € [—-x/3,x/3]. 


Au regard des spécificités des caméras embarquées sur les drones [56], nous constatons que le 
positionnement de celles-ci a peu d’influence sur la trajectoire du vecteur. 

Dans les simulations correspondant aux Figures 3.4 et 3.5, nous avons choisi des caméras de ce 
type. 
Remarque 3.3. Dans la suite du travail, nous considérons uniquement le problème de planification 
de trajectoires du vecteur. En effet, comme la caméra embarquée peut suivre la cible sans difficulté, 
nous ne traitons que la planification du vecteur car cela est suffisant pour éviter que la cible ne sorte 
du champ de vision de la caméra (e.g., à cause de la limite de zoom, de l’occlusion de la cible...). 


3.1.3 Découplage de la caméra 


Dans le paragraphe précédent, nous avons conclu que la problématique du couplage vecteur/capteur 
se résume au guidage du vecteur (les caméras utilisées étant très performantes et non contraintes). 
C’est pourquoi, dans cette partie, nous proposons une méthode de commande indépendante de la 
caméra, que le problème soit 2D ou 3D. 


3.1.3.1 Étude du cas plan 


Nous modélisons le vecteur par le système de Dubins (2.2) et, comme au paragraphe 3.1.1, nous 
supposons que la caméra agit comme un système linéaire du premier ordre et qu’elle est commandée 
par une commande %e, i.e. 


1 
& = u + —uc — —(a — 6) 
Te Te 


avec «à définissant l’angle de la caméra dans le repère global, i.e. à = & + 0 où à est l’angle de la 
caméra dans le repère associé au drone. La caméra est supposée non contrainte, i.e. à € R. 
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CHAPITRE 3. RÉSUMÉS DES CONTRIBUTIONS PERSONNELLES 


3.1. Comment effectuer le couplage vecteur/capteur ? 


ñ 1 ñ ñ 4 1 ni ñ 1 1 n 1 1 1 
0 5 10 15 20 0 5 10 15 20 30 35 40 45 50 


25 
x Temps 
(a) |[El| = 1:3091, oyey = 0.6500 (b) & = —0.6579, 0. = 1.1576 


FIGURE 3.3 — Trajectoire primant l’erreur de positionnement du vecteur. Le champ de vision de la 
caméra est supposé contraint : à € [—-7/3,7/3]. 


-0.05 F 


TS 
S 015 
© 
w -0.2 
el -0.25t 
A4f L L L L | L L L L L L L 
0 5 10 15 20 25 03 5 10 15 20 25 30 35 40 45 50 
x Temps 
(a) [ll = 1.2232, ojey = 0.6040 (b) & = —2.9664e—4, ©. = 0.0093 


FIGURE 3.4 — Trajectoire primant l'erreur de positionnement de la caméra. Le champ de vision de 
la caméra n’est pas contraint : & € R. 
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CHAPITRE 3. RÉSUMÉS DES CONTRIBUTIONS PERSONNELLES 


3.1. Comment effectuer le couplage vecteur/capteur ? 


-0.05 F 


>, Ca) 
S 015 
al ue” 
Ô 
al © W 52! 
0 
0.25} 
_2 
0 5 10 15 20 086 5 10 15 20 25 30 35 20 45 50 
x Temps 
(a) |[El| = 1:3092, ojey = 0.6507 (b) & = —2.9495e—4, 0e, = 0.0093 


FIGURE 3.5 — Trajectoire primant l’erreur de positionnement du vecteur. Le champ de vision de la 
caméra n’est pas contraint. 


JHANAnITINAN 
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| PUUUVUT Ty 
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(b) & = —0.1263, 0e, = 0.7688 
FIGURE 3.6 — Sur cette figure les mouvements de la caméra sont contraints : à € [-7/3,7/3]. 


L'objectif est que la caméra suive le pattern en trait pointillés tout en assurant que la trajectoire ne 
franchisse pas la droite verticale. 
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CHAPITRE 3. RÉSUMÉS DES CONTRIBUTIONS PERSONNELLES 


3.1. Comment effectuer le couplage vecteur/capteur ? 


Nous supposons connues la position et la vitesse de la cible, respectivement qcible €t dipl. Comme 
défini précédemment, nous souhaitons asservir l’angle de la caméra sur celui nous permettant de 
suivre (visuellement) la cible choisie. Notons @nr Cet angle, nous avons : 


ju Ucible — Y 
Qcible = AICtan | —————— 
cible — 


ag, = Cible  )(rcibie — T) — (cible — Y)(Hciie — À) 
D (Tcible — T)? + (Ycible — Y)? 


Remarque 3.4. Si tabie — æ = 0 alors nous posons @cible = 7. Si (Zciple — €)? + (Ycible — Y)? = 0, nous 
considérons cible = 0. 


Étant donné que l’équation de la caméra est un système linéaire d’ordre 1, nous considérons que 
l'erreur décrit un système différentiel d’ordre 1, i.e. l'erreur &; = @cipie — @ vérifie 


Êe = LEe, 


avec & < Ü une constante réelle qui nous permet d'obtenir la stabilité de la solution stationnaire 
Ee = 0. 
À cet effet, nous proposons la commande suivante [102] 


| 1 
Uc — Te (ane = ut — (a — 0) — se) ; 


C 


3.1.3.2 Extension au cas 3D 


Le cas d’un système évoluant dans l’espace à trois dimensions est légèrement différent, mais le 
principe reste le même. 

Considérons le système d’équation (2.6), modélisant le drone MALE, comme présenté au para- 
graphe 2.1.2. 

Nous supposons que la caméra à deux degrés de liberté en rotation qui nous permettent de 
modifier l’azimut 4 (la rotation autour de son axe vertical, parallèle à l’axe k) et l’angle d’élévation 
de la direction d'observation 4, voir Figure 3.7. 

Comme au paragraphe 3.1.1, nous supposons que la caméra est modélisé par deux systèmes 
linéaires du premier ordre découplés (de constantes de temps 7, et 7. et de gains unitaires) et 
commandés par up, et uy,, 1e. 


- 1 1 
be = —u9, — —06. 
TB. TO. 
- 1 1 
Ve = pe — —%Ÿ, 
The The 


Nous souhaitons asservir la direction d'observation de la caméra à celle nous permettant de 
suivre (visuellement) la cible choisie. La vitesse et la position de cette dernière sont supposées 


50 


CHAPITRE 3. RÉSUMÉS DES CONTRIBUTIONS PERSONNELLES 


3.1. Comment effectuer le couplage vecteur/capteur ? 


m7" | 


(4 i 


FIGURE 3.7 — Schéma du drone et de sa caméra en 3D 


connues. Cette direction à une composante d’azimut 4e et une d’élévation Vin. Nous obtenons 
ces composantes ainsi que leurs dérivées de la manière suivante, 


AVE 
Ocibre = arctan | ——= 
cible (2) 


à : = Ax2AT: Fr AxtoAT: 
cible — (Az)? Fe (Az)? 


et 


VYible = arctan — 
(Ari)? + (Ax2)? 


((Ax:)? + (Azx2)°) ES — AÂx3 (ann: arr) 
PR rm a) ae Ga) 


avec At; = Tip — Li Et Li, Tiavrer à = {1,2,3} les coordonnées, respectivement, de la caméra et de 
la cible dans le repère global. 


Remarque 3.5. Si Az; = 0 alors nous posons cible = T. Si (Az)? + (A2)? = 0, nous considérons 
VYcible = ñ/2, VYcible =0et Ocibie = 0. 


La caméra étant un système linéaire d'ordre 1, nous considérons que l'erreur décrit aussi un 
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CHAPITRE 3. RÉSUMÉS DES CONTRIBUTIONS PERSONNELLES 


3.2. Contrôle de type Lyapunov et temps-minimal pour les drones 


système différentiel d'ordre 1, i.e. les erreurs £g, = bcibie — 0e et Em. = Vcible — Ye Vérifient 


ÉD. — LED, 
Eve — LEe 


avec & < 0 une constante réelle qui nous permet d’avoir la stabilité des solutions désirées. 
Ces objectifs sont réalisés en utilisant les commandes suivantes [102] 


. 1 
UD. = TE. (é + — 04 — LEg. 
TO 


(a 


; 1 
Uype — The (an + a _ as) 


3.2 Contrôle de type Lyapunov et temps-minimal pour les drones 


L'article présenté dans ce paragraphe est en cours de publication. 
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Lyapunov and minimum-time path planning for drones 


Thibault Maillot* Ugo Boscainii Jean-Paul Gauthier“ Ulysse Serresÿ 


July 18, 2013 


Abstract 


In this paper, we study the problem of controlling an unmanned aerial vehicle (UAV) 
to provide a target supervision and/or to provide convoy protection to ground vehicles. 

We first present a control strategy based upon a Lyapunov-LaSalle stabilization method 
to provide supervision of a stationary target. The UAV is expected to join a pre-designed 
admissible circular trajectory around the target which is itself a fixed point in the space. 

Our strategy is presented for both HALE (High Altitude Long Endurance) and MALE 
(Medium Altitude Long Endurance) types UAVSs. A UAV flying at a constant altitude 
(HALE type) is modeled as à Dubins vehicle (i.e. a planar vehicle with constrained turning 
radius and constant forward velocity). For a UAV that might change its altitude (MALE 
type), we use the general kinematic model of a rigid body evolving in R°. Both control 
strategies presented are smooth and unlike what is usually proposed in the literature these 
strategies asymptotically track a circular trajectory of exact minimum turning radius. 

We also present the time-optimal control synthesis for tracking a circle by a Dubins 
vehicle. This optimal strategy, although much simpler than the point-to-point time-optimal 
strategy obtained by P. Souères and J.-P. Laumond in the 1990s (see [45]), is very rich. 

Finally, we propose control strategies to provide supervision of a moving target, that 
are based upon the previous ones. 


Keywords: optimal control, path planning, aircraft navigation, unmanned aerial vehicles, 
rigid-body dynamics, under-actuated systems, nonlinear control, trajectory tracking. 


1 Introduction 


This study holds in the context of the French FUI SHARE project (see [3]), supported by a 
consortium of companies and research labs !, and the authors are granted from the project. 
The aim is to design a trajectory for a UAV from some starting point and orientation to some 
final submanifold of the state space. 

The UAV path planning problem is a standard motion planning problem and although the 
bibliography on this subject is very rich we do not aim to provide an exhaustive list. Let 
us just mention that classical planning algorithms are based upon very different approaches 
such as geometric control [14, 19, 32], optimal control [21], flatness [7], stochastic theory [5]. 
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Saclay, 91128 Palaiseau Cedex, France; email: ugo.boscain@polytechnique .edu 
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SUniversité de Lyon, F-69622 Lyon, France; Université Lyon 1, Villeurbanne; LAGEP, UMR CNRS 5007, 
43 bd du 11 novembre 1918, 69100 Villeurbanne, France; email: ulysse.serres@univ-lyoni.fr 

1 Opéra Ergonomie, ONERA, Thales Alénia Space, Eurocopter, Adetel group. 


Moreover, the reader may refer to [27] for a survey on trajectory-planning algorithms focusing 
on autonomous UAV guidance. 


The purpose of this paper is to present different solutions to the above mentioned motion 
planning problem for both HALE and MALE types of fixed-wing UAVSs. Although the case 
of rotary-wing drones is interesting, it is simpler and more standard motion planning methods 
may be applied (see [8, 26]). 


The problem is addressed from a kinematic point of view only. In particular, we consider 
that the aircraft controls are just the angular velocities. From a kinematic point of view, a 
UAV can be specified by a control system q = f(q,u) (specific models will be given later on), 
where q denotes the state of the UAV which usually includes the UAV position (in R? or R°) 
and attitude (ie. the UAV orientation in S1 or SO(3)), and u is the control vector driving 
the UAV kinematics, such as angular and spacial velocities. Moreover, we make the following 
assumptions on the UAV: 


e the velocity of the UAV is assumed to have à positive lower bound (no stationary or 
quasistationary flights are allowed) and a positive upper bound; 


e the UAV is assumed to be kinematically restricted by its minimum turning radius r > 0, 
or equivalently, its yaw angle is assumed to be constrained by an upper positive bound. 


In this article, tools from Lyapunov-based stability analysis and optimal control theory are 
used for the development of control algorithms achieving the desired motion planning task. 

The paper is organized as follows: in Section 2, we present global Lyapunov-LaSalle-based 
stabilization results for a UAV tracking a horizontal circle of minimal turning radius. Section 
2.1 and Section 2.2 are the detailed studies of the cases of HALE and MALE UAVSs respectively. 

Section 3 describes the time-optimal synthesis for the Dubins system tracking a minimal 
radius circle. The study is done in a reduced state space. On this subject, we have to mention 
the great work [45] about the more complicated point-to-point problem (see also [1] for a partial 
study). 

Finally, in Section 4 the results of Sections 2 and 2.2 are brought together and an algorithm 
for the tracking of a moving target is proposed. In the present paper, what is meant by target 
is the position (in R? or RŸ) of the center of mass of the object the UAV has to follow. 


We end this introductory section providing a list of mathematical conventions used here in. 
In the paper, smooth means C'*. Throughout the remaining of the paper, the time dependence 
of variables is often omitted to lighten notations. Moreover, we use 


e the dot ‘ to denote the derivative with respect to time; 


[x] to denote the usual Euclidean norm of x € R”, and (x, y) for the inner product of two 
such vectors; 


B(x,r) € R” to denote the open ball of radius r centered at x; 
e M! to denote the transpose of a matrix M; 
e tr(M) =), M; to denote the trace of a n x n square matrix M; 


e [:,:] to denote the Lie bracket between vector fields and the commutator of matrices. 
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e For two n x n real matrices M and N, the Hilbert-Schmidt scalar product is tr(M'N), 
then the Hilbert-Schmidt norm of a matrix M is /tr(M'M). 


e Also, we do not distinguish between row and column vectors, the distinction being clear 
from the context. 


2 Lyapunov-LaSalle-based stabilization 


This section is concerned with the problem of global stabilization of a UAV (in the position- 
attitude state space) on à fixed horizontal circle of minimum turning radius with prescribed 
direction that we shall refer to as the final manifold. 

By using a standard Lyapunov-LaSalle-based approach, we propose, for both HALE and 
MALE type UAVs, new nonlinear feedback control laws that stabilize the UAV on the final 
manifold. 

The Lyapunov method has been used by several authors, in different ways such as for 
example [15, 16, 20, 22, 24, 28, 33, 35, 50]. Here we present an original Lyapunov strategy, 
serving the advantage to be completely smooth, very simple to apply, and with the peculiarity 
that the circular final manifold is with exactly minimum turning radius. Frequently, Lyapunov 
strategies presented in the literature do not have this last feature. 

Along this section, controls are assumed to be smooth. Nonsmooth controls appear in 
Section 3 only. 


2.1 Lyapunov based method for the Dubins system 


2.1.1 Problem under consideration 


From the kinematic point of view, a rough HALE drone is governed by the standard Dubins 
equations (see e.g. [18] for a justification, for Dubins problems on Riemannian manifolds see 


[17, 43]): 


à = cos Ü 
ÿ = sin 0 (2.1) 
Ô = u. 


with (x,y,0) € R? x S! being the state (where (x,y) € R? is the UAV’s coordinate in the 
constant altitude plane, and 0 the yaw angle), and u € [— max, Uumax] being the control variable. 
Note that the yaw angle 0 is the angle made by the aircraft direction with respect to the x-axis. 

These equations express that the drone evolves on a perfect plane (perfect constant altitude), 
at perfect constant speed 1, moves in the direction of its velocity vector, and is able to turn 
right and left. 

Note that the rough model (2.1) is pertinent for HALE drones, the velocity being almost 
really constant. 

As said in the introduction, the UAV is assumed to be kinematically restricted by its positive 
minimum turning radius r which implies the following bound on the yaw angular velocity: 

il 
Tr = max SU < Umax = 

Remark 2.1. Note that up to a dilation in the (x,y)-plane we may assume without loss of 
generality that [—unax, Umax] = [-—-1,1]. This normalization could be used to simplify the 
treatment. We do this in Section 3. 
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We denote by C the final manifold which is defined to be the counterclockwise oriented circle 
centered at the origin of the (x,y)-plane, with radius r = 1/umax. In the (x,y,0)-coordinates 
system, C is given by 

C={(x,y,0)| x =rsin0, y = —-rcosb}. 


We have the following result. 


Theorem 2.2. There exists a smooth feedback control k : R? x S! — [—-Umax; Umax] SUCh that 
the set C is a globally asymptotically stable attractor for the closed-loop system resulting from 
applying the feedback control u = k(x,y,0) to system (2.1). 


2.1.2 Proof of theorem 2.2 


The proof of Theorem 2.2 is based upon a standard Lyapunov-based control design and LaSalle’s 
principle (see e.g. [25, Theorem 3.5, p.190] or the original LaSalle’s paper [31]). 

To construct a state feedback stabilizer to € for system (2.1) (and further to construct the 
time-optimal synthesis in Section 3), it is convenient to use UAV-based coordinates (£, ÿ, 0) 


with % and ÿ defined by 
&\ cos 0 sin 0 zx 
ÿg)/  \-sin0 cos0) \y}° 


in which the set C takes the form 
Ê= {(&, 9,0) | T = 0, y = m. 


Notice that if C would have been travelled clockwise, the equation ÿ = —r would have been 
changed for ÿ = r. 

In the (#,ÿ,0)-coordinates system, the variable 0 is decoupled and the two first equations 
of system (2.1) may be rewritten as 


À = uÿ +1 
Re (2.2) 
ÿ = —u. 


Remark 2.3. One can read, in (2.2) that if u is changed for —u, the two equilibria of (2.2) 
corresponding to the control values —umax and Umax are exchanged so that the set of equilibria 
remains unchanged. It means that we can indifferently consider one of the two equilibria 
positions. 


Changing the (#,ÿ) coordinates for (&,ÿ) = (&,ÿ + r), the stabilization problem to C is 
reformulated in these variables as the stabilization problem to the submanifold {(%, ÿ,0) | x = 
ÿ = 0}. Equivalently, we will refer to the convergence of (%, ÿ) to the point (0,0) of the reduced 
state space. The reader can easily check that the (%,ÿ) obey the following equations: 


(2.3) 


The proof of Theorem 2.2 is a consequence of its analogue in the reduced (%, ÿ)-coordinate state 
space. 


Lemma 2.4. There exists a (smooth) feedback control k : R? — [-Umax: Umax] which stabilizes 
trajectories of the system (2.3) with u(t) = k(x(t), y(t)) with respect to the origin. 
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Proof. We shall prove that 
VGD=r+T 

is à Lyapunov function for system (2.3) if the control law is well chosen. A straightforward 
computation yields 

V(&,ÿ) = 27 (1- ru). 
To be a Lyapunov function for system (2.3), we need to ensure that V is nonpositive. Keeping 
in mind that Umax = 1/7, it vields u = umax when æ > 0. Let € be a positive real number, let 
a be such that —Umax  @ < Umax and let & : R? — [-Umax; Umax] be à (actually any) smooth 
feedback control which satisfies 


RG ) Umax fxz>0 (24) 
3,ÿ) = , . 
: —Umax < k(&, y) << —Umax fx <0. 


For instance, one can take the function defined by 


a if x < —€ 
k(t,y) = anse +0 ifxe(—-Ee,0) 
Umax fx>0, 


the shape of which is drawn on the next figure. 


k (x) 


Umax 


According to LaSalle’s principle and since V is a proper function, we infer that all the trajec- 
tories of system (2.3) with feedback control k(-) converge to the largest invariant set contained 
in the set E = {(%,ÿ) | V(&,ÿ) = 0} = {(%,ÿ) | 3 > 0}. Due to condition (2.4) on the feedback 
control, it is easy to see that any solution to (2.3) starting at (Zo, Yo) with Zo > 0 escapes from 
E in finite time. Indeed, when xo > 0 the solution to the closed-loop system satisfies on E: 


= UmaxŸ 
ÿ = —UmaxT. 
Hence, the solution (x(t),y(t)) = (to CoS(ümaxt)+ÿ0 SIN(Umaxt), Yo COS(Umaxt))— To Sin(Umaxt) 
escapes from Æ at time 


1 5 __ 

a arctan (2) if Yo < O 

t — = if ÿo = 0 
1 à __— 

 . fr — arctan (æ)) if Yo > 0. 


Gi 
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Consequently, the largest invariant set contained in Æ is the equilibrium point (0,0). This ends 
the proof. HE 


Proof of Theorem 2.2. Since V(x,ÿ) tends to zero implies that (x,y,0) tends to C, Lemma 
(2.4) shows that k steers asymptotically system (2.3) to C, or equivalently that k(x,y,0) = 
k(o(x, y, 0), (where ® : (x,y,0) + (&,ÿ)) steers asymptotically system (2.1) to C. = 


2.1.3 Numerical results 


The theory exposed in the preceding subsections was tested by means of simulated data. During 
these simulations, we numerically integrate the closed-loop system resulting from applying the 
feedback control u = k(x, y, 0) of Theorem 2.2 to system (2.1) during 200 units of time. 

Figure 1 shows two trajectories of the same UAV starting from the same initial point but 
with two different values of the parameter € in the control function defined by equation (2.4). 
On both pictures the UAV reaches the final manifold tangentially, here the circle of radius 2 
centered at the origin. 


ec 


Figure 1: Illustration of Lyapunov-LaSalle’s-based 2D-planification results with starting point 
(to, Yo, 00) = (20, 0,0), control constraints max = 0.5 and a = 0.1. (a) € = 2, (b) € = 0.5. 


2.2 Extension of the method to the 3D-case 
2.2.1 Statement of the problem and of the result 


In this section, we extend the results of the previous section to the case of MALE UAVSs flying 
at constant speed. Our goal is unchanged: we want to (globally) stabilize the UAV on a final 
manifold. 

From the kinematic point of view, a rough MALE drone (modeled by a standard 6 degree of 
freedom solid, symmetric w.r.t. its vertical middle plane) is governed by the following equations 


(2.5) 


X = RV 
R= RQ 


with 
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e Q=(X,R) € R° x SO(3) being the state where X is the UAV's coordinates expressed 
in the inertial frame and R, the UAV attitude, is the rotation matrix between the UAV’s 
fixed frame and the inertial frame which gives the motion of the UAV’s frame orientation; 


e V5 = (1,0,0) being the thrust direction expressed in the UAV's frame; 


O0 —u3 
e and ( — ( us 0 #) € 50(3) (the Lie algebra of SO(3) consisting of skew-symmetric 


—U2 U] 0 
matrices) the angular velocity of the body expressed in the body-fixed frame. 


System (2.5) is used by many authors in the literature to model UAVS and airplanes (see e.g. 
[29, 48]. See also [49] and references therein. We define u = (u1,u2,u3) to be the control 
variables and we denote by U — EME Ui max] the control set. The control variables are 
the roll (uw), pitch (u2) and yaw (u3) angular velocities and are constrained by ü;max > 0 for 
i = 1,2,3. 


Remark 2.5. Notice that the controllability of system (2.5) with control in U is à direct 
consequence of [9, Theorem 1] (see also [40, Theorem 6.7]). 


Equations (2.5) express that the drone moves in a 3-dimensional space, at perfect constant 
speed 1, moves in the direction of its velocity vector, and is able to turn in every direction of 
its fixed frame. 

The First equation of system (2.5), in R°, describes the evolution of the UAV in space. 
The second one, in SO(3), depicts the motion of the UAV’s frame which is controlled by three 
controls on angular velocity acting physically on the roll, pitch and yaw angles. 

Set X = (æ1,%2,%3). The (x1,%2)-plane is called the horizontal plane. Here, we fix the 
final manifold C to be the union € = C1 UC2 of the two circles of radius r3 = 1/u3max Centered 
at the origin in the horizontal (x1,æ2)-plane and oriented counterclockwise in the orientation 
given by their normals M = (0,0,1), N2 = (0,0, —1). Define 


00 0 0 O1 0-10 
J=1|00-1|, 2=| 000, Z=1|11 0 0), 
01 0 —100 0 0 0 
so that Q = u1J1 + u2Jo + u3J3. 
Notice that the Lie algebra structure on 50(3) is given by 
(a, do] = ds, [o,J3]= Ji, [s3, di] = da. (2.6) 


The reader can easily check that, in the (X, R)-coordinates system, the final manifold takes 
1 0 0 
the form: € = C1 U Co, setting A; = (o 6). 


C1 = {(X,R) € R° x SO(3) | X = r(sin0,—cos0,0), R=e/#, 6ER} 


C2 = {(X,R) € R° x SO(3) | X = ra(sin0,cos 0,0), R = Aje/, O6ER}. 
We will prove the following result. 
Theorem 2.6. There exists a (smooth) feedback control k : R° x SO(3) — U such that the set 
C is a global asymptotic attractor for the closed-loop system resulting from applying the feedback 
control u = k(X,R) to system (2.5). Moreover, C1 is an asymptotically stable limit cycle and 


C2 is an unstable limit cycle. AU trajectories tend to C1 but the trajectories starting in the plane 
x3 = 0 with R(0) = Ae”30 for some 0, which tend to Co. 
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2.2.2 Proof of Theorem 2.6 


Here again, to construct a state feedback stabilizer to C for system (2.5) it is convenient to use 
UAV-based coordinates (X, R) defined by X = RX in which X obeys 


> ee 
and the final manifold takes the form € = a U Ca; 
G = {(X,R) € R* x SO(3)| À = (0,-rs,0), R=e*?, OER}, 


= {(A,R) € R° x SO(3) | À = (0,-r3,0), R= Aie, 0€ R}. 


Set X* = (0, —r3,0) and consider the new variable 


It vields the new system 
(2.7) 


In particular, we have 
Zi — U3T2 = Uu2T3 FF (1 — Tau3) 


To = —U3%1 + T3 


ä3 = UT) — U1T2 + Tau]. 


The stabilization problem to C is reformulated in these variables as the stabilization problem 
to the submanifold € = C1 UC», 


Ci ={(X,R) ER° x SO(3)| ZX =0, R=e"#, 8ER}. 


Co={(X,R)ER° x SO(3)|X =0, R=Ae/#, BER}. 


Let o be a positive real number. We shall prove that the function V : R° x SO(3) — R, 
defined by 


VIE, R) = : (oX'X + tx ([R, JR, Ja) 
= Vi(X, R) + V(X,R), 


with Vi(X,R) = SX'X and Vo(X,R) = itr([R, B3J'[R, Ja], is a Lyapunov function for the 
closed-loop system resulting from applying a well-chosen feedback control u = k(X, R) to (2.7). 


(2.8) 


Remark 2.7. Notice that the coefficient o in the Lyapunov function (2.8) is a weight parameter 
that can be tuned in order to give more (resp. less) importance to the horizontality (i.e., pitch 
and roll angles close to zero) of the drone rather than its distance (in R°) to the target if o < 1 
(resp. o > 1) during the stabilization process. Although this parameter plays no role in the 
subsequent results and computations, we keep it in our formulas. 


The function V is clearly smooth. Moreover, the following proposition is obvious. 


Proposition 2.8. V(-) is a proper function. 
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Notice that V (X, R) may be rewritten as 


V(X,R) = XX +2+tr(RJRJ3). 


WI a 


Hence, taking into account equations (2.6), the derivation w.r.t. time of V (X,R) along the 
trajectories of system (2.7) vields 


V (Z,R) = oX'X + Sur (RJ3R'J3) 


oX'(Vo — Q(X + X*)) +tr (R[Q, Js]R'J3) 
= oX' (V6 - QX*) + tr(R(ui[J1, Ja] + uo[Jo, J3]R)J3) 


T1 1 0 

© a |, [0] +(uidi +u3ds) | ra ) — ujtr (RIR!J3) + uotr (RJR!J3) 
T3 0 0 

= où (1— raus) + (ora%3 — tr (RJ2R'J3)) + uotr (RJ1R'J3) . 


Remark 2.9. Note that Vi (X,R) = o%1 (1— raus) and Vo (X,R) = à (orats — tr (RJRJ3))+ 
uatr (RJ1R'J3). Later it will be important that both quantities will be negative. 


Let € be a (small) positive real number. Set £ = orst3-tr (RBR'J3) and & = tr (RJ1R'J3). 
In order to ensure the non positiveness of V, we choose smooth feedback controls u1 = k1(X,R), 
u2 = k(X, R) and u3 = k3(X, R) such that: 


sign k1(X, R) = — sign £ (2.9) 
sign k2(X, R) = — sign (2.10) 


= 3 max if 7 2 0 
ka(X, R) = { Rue (2.11) 


—U3 max < k3(X, R) U3 max if Zi < 0. 
The feedback controls k1, k2 may be rewritten as 


k1(X,R) = —E1çg1 (1), 
ko(X, R) = —£2po (é2), 


where @1, 2 are smooth strictly positive functions and £1 (resp. £2) is zero if and only if u, = 0 
(resp. u2 = 0). Hence, 


V(X,R) = o%1 (1 — raus) — Éfp1 (£1) — 40 (4). 


According to LaSalle’s principle and since V is a proper function, we infer that all the trajec- 
tories of system (2.7) with feedback controls k1, k2 and k23 converge to the largest invariant set 
T contained in the set 


E = {(X,R)| V(X,R) = 0} 
un {(X,R) | k(X,R) = ko(X, R) = 0, k3(X, R) = U3 max } 
= {(X,R) | êi = £2 = 0, k3(X, R) . U3 max }- (2.12) 


Obviously we have € € I. Let us prove the converse. Due to conditions (2.12) on the 
feedback controls, it is easy to see that any solution to (2.5) starting from (X,R) € E and 
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remaining in Æ satisfies 
T1 = U3 max T2 
L2 = —U3maxT1 
ns (2.13) 
T3 — 0 
R = Rus maxJ3. 


Consequently, any solution to (2.13) starting from 7 with %1 > 0 escapes from E in finite time. 
Thus, any point in 1 satisfies T1 — 0 and consequently, due to the first equation of system 
(2.13), ï2 = 0. Moreover, according to the third equation of system (2.13) Z3 is constant in 
E. This constant is easily seen to be equal to zero. Indeed, let t++ (X(t), R(t)) be a smooth 
curve in Î defined on an interval around zero. Since the functions £1 and &2 are identically zero 
in E, we have 


OT3T3 = tr (RIR'J3) : (2.14) 
O=tr(RJR'J3). (2.15) 


Since #3 is constant, differentiation of equation (2.15) (divided by ü3max) W.r.t. t yields 


“50 
0 = er (RJR J3) 


— : tr (Rlu3 max J3, Ji]R'J3) 


U3 max 


= tr (RJ2R'J3) , (2.16) 


which, according to (2.14), implies that Z3 = 0. Summing up, we get 1 € {71 = Z2 = 3 = 0}. 
It remains to show that À has the required form. Equations 2.15 and 2.16 may be rewritten 
as tr(Zi(R'J3R)) = 0 and tr (J2(R'J3R)) = 0 showing that R'JR is orthogonal to Ji and Ja 
for the Hilbert-Schmidt scalar product. It follows that R/J3R = \J3 with À € R. Moreover, 
conjugations by elements of SO(3) being isometries for the Hilbert-Schmidt scalar product, we 
must have |A] = 1. Consequently, R'J3R = +J3. The case À = +1 exactly means that the 


rotation matrix À commutes with J3. According to Lemma 2.10 below, R is of the form es. 
J30 


The case À = —1 is similar and we conclude that R is of the form Aie 
It is clear that all trajectories 7 (t) starting in the plane {x3 = 0} with R(0) = Ae”#? tend 
to C2 since uw = 0, u2 = 0. 
Let us show that any other trajectory tends to C1. Due to Remark 2.9, V> < 0 then W 
cannot increase along any trajectory. But V2 reaches its maximum exactly on the union of the 
trajectories of type T(t). 


Lemma 2.10. The stabilizer S3 of J3 under conjugation by SO(3) is the subgroup of rotations 
ï 0 sin0 0 
R3 about the x3-axis, R3 = (= cosb 6). 


Proof. With obvious block notations, write J3 = (49) and H = (4). Then, [J3,H] = 
(es nr For H € S3, H'J3H = J3, which is equivalent to [J3, H] = 0. It implies that 


B =0,C =0 and [J, A] = 0. Since H € SO(3), we first conclude that d = +1 and À € O(2). 
The stabilizer of J under conjugation by O(2) is SO(2). Then, À € SO(2), det(A) = 1, 
[ 


and d = 1. Therefore, H is of the form R3. 
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2.2.3 Numerical results 


The theory exposed in the preceding subsections was tested by means of simulated data. During 
these simulations, we numerically integrate the closed-loop system resulting from applying the 
feedback control u = (ui,u2,u3) = (k1(X,R),k2(X, R), k3(X, R)) of Theorem 2.6 to system 
(2.5) during 200 units of time. 

The feedback control k3 is built as its analogue k in the 2D-case. Concerning the feedback 
controls k1 and k> we chose smooth functions varying slowly between the minimum and the 
maximum control constraints. 

Figures 2 and 3 display trajectories obtained for the same UAV which starts at the point 
(50,20,20) with a direction of 0 radian for the pitch, roll and yaw angles. We consider the 
following curvature constraints, Ujiimax = U2max = —Uimin = —U2min — 0.1, U3max — 0.5, 
U3 min — 0.1. 

The difference between Figure 2 and Figure 3 is due to the value of the parameter € in the 
control function defined in equation (2.11). In the first case we set € = 3 and we use the value 
€ = l'in the second simulation. In each case we set & = 1 in (2.8). 


Figure 2: Illustration of Lyapunov-LaSalle-based 3D-planification results with starting point 
(t10, 20, Z30) = (50,20, 20), control constraints U1 max = U2max — 0.1, Ugÿmax = 0.5, az = 0.1, 
o = lande =3. (a) The (æ1(-),x2(-),x3(-)) trajectory, (b) projection in the horizontal plane 
(left) and evolution of z3 w.r.t. time (right). 


3 The time-optimal stabilizing synthesis 


In this section we study the minimum time problem for system 2.1. For simplicity of exposition 
we assume that [—Umax, Uümax] = [—1,1] (cf. Remark 2.1). More precisely, we consider the 
following problem that we shall denote by (P). 


(P) For every (0, Yo, 00) € R? x S! find the pair trajectory-control joining (xo, Yo, 60) to €, 
which is time-optimal for the control system 


à = cos Ü 
ÿ=sin@ (3.1) 
Ü=u, ue [1,1]. 


WT 
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0 10 20 30 40 50 
æi 


(a) (b) Projection in the horizontal plan and evolution of 
æ3 over the time. 


Figure 3: Illustration of Lyapunov-LaSalle-based 3D-planification results with o = 1 ande = 1. 


To solve Problem (P) it is again convenient to work with a reduced system in dimension two. 
Indeed, in dimension two, a complete theory for finding the time-optimal synthesis exists and 
will be recalled in the next two sections. In this section we make the following transformation 


(in O(2) \ S0(2)): | _—. 
(5) = Ge ee) 6) 


This transformation decouples the variable 0 and the final manifold C projects to the point 
(&,ÿ) = (0,1). We thus consider the following problem that we shall denote by (P’). 


(P’) For every (Zo,ÿo) € R? find the pair trajectory-control joining (Zo, go) to go = (0,1), 
which is time-optimal for the control system 


The collection of all solutions to Problem (P’) for every (o, ÿo) is called the time-optimal 
stabilizing synthesis ([371). 

First, it is convenient to change the sign of the dynamics and to consider the equivalent 
problem (denoted by (Q)) of finding the time-optimal synthesis issued from qo: 


(Q) For every (%,ÿ) € R? find the pair trajectory-control joining go = (0,1) to (&y,ÿr), 
which is time-optimal for the control system 


Fe ue [1,1]. (3.2) 


Once that Problem (Q) is solved, then the time-optimal stabilizing synthesis is obtained simply 
by inverting the arrows of the trajectories of the time-optimal synthesis. 
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3.1 Basic definitions and Pontryagin Maximum Principle 


Let M be a n-dimensional connected smooth manifold and let F and G be smooth vector fields 
on M. Consider the following general control-affine time-optimal problem that we shall denote 


by (R). 


(R) For every qo and qf in M find the pair trajectory-control joining qo to gr, which is time- 
optimal for the control system 


g=F(g)+uG(qg), qe M, uel-1,1l]. (3.3) 


Definition 3.1 (admissible control/trajectory). An admissible control u(-) for the system (3.3) 
is a measurable function u(-) : [a,b] — [—-1,1]. An admissible trajectory is a Lipschitz function 
g(-) : [ab] — M satisfying q(t) = F(q(t)) + u(t)G(q(t)) a.e. for some admissible control u(-). 


In the following we assume that the control system is complete ï.e., for every measurable 
control function u(-) : [a,b] — [1,1] and every initial state go, there exists à trajectory q(:) 
corresponding to u(-), which is defined on the whole interval [a, b] and satisfies q(a) = go. 

Thanks to the compactness of the set of controls, the convexity of the set of velocities and 
the completeness of the control system, Filippov’s theorem (see for instance [1]) gives, 


Proposition 3.2. For any pair of points in M, there exists a time-optimal trajectory joining 
them. 


The main tool to compute time-optimal trajectories is the Pontryagin’s Maximum Principle 
(PMP for short). We refer to [1] for a general version of PMP and also [12] for a version of 
PMP on two-dimensional manifolds for single input control affine systems. The PMP is a first 
order necessary condition for optimality. 

The following theorem is a version of PMP for control systems of the form (3.3) that we 
state in our own context only. 


Theorem 3.1 (PMP for problem (R)). Consider the control system (3.3). Define for every 
(g,p,u) € T*M x [-1,1] the function 


H(q,p,u) = (p, F(q)) + u (p,G(q)). (3.4) 


If the pair (q(-),u(-)) : [0,T] — M x [-1,1] is time-optimal then there exists a never-vanishing 
Lipschitz covector p(-) : t € [0,T] + p(t) TM and a constant À < 0 such that for a.e. 
te [0,T]: 
0H 
. 400) = (DD (D), 


ii. p(t) = CO OO) 


ii. H(q(t),p(t),u(t)) = en H(q,p,u), 


iv. H(q(t),p(t),u(t)) + À = 0. 


Remark 3.3. À trajectory q(-) (resp. a pair (q(-),p(-))) satisfying the conditions given by the 
PMP is said to be an extremal (resp. an extremal pair). An extremal corresponding to À = 0 
is said to be an abnormal extremal, otherwise we call it a normal extremal. 
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3.2 An overview on optimal synthesis on 2D-manifolds 


In this section we briefly recall the theory of time-optimal syntheses on 2D-manifolds for systems 
of the form (3.3), developed by Sussmann, Bressan, Piccoli and Boscain in [11, 13, 36, 46] and 
rewritten in [12]. 

This section is written to be as much self-consistent as possible. In the remainder of the 
section, we assume M to be of dimension two. 


3.2.1 Basic definitions 


For every coordinate chart on the manifold M it is possible to introduce the following three 
functions: 


Aa(g) = det(F(q), G(g)), (3.5) 
Ag(g) = det(G(g),[F, G1(g)), (3.6) 

and on M \ A,1(0) r 
fs(o = - RU (37) 


The sets A31(0), A5! (0) of zeros of A4, Ag are respectively the set of points where F and 
G are parallel, and the set of points where G is parallel to [F, GT. These loci are fundamental 
in the construction of the optimal synthesis. In fact, following [12], assuming that they are 
smooth embedded one dimensional submanifolds of M we have the following: 


e on each connected component of M \ (A3!(0)[UJA%'(0)), every extremal trajectory is 
bang-bang with at most one switching. Moreover, for every switching of the extremal 
trajectory the value of the control switches from —1 to +1 if fs > 0 and from +1 to —1 
if fs < 0; 


e the support of singular trajectories (i.e., trajectories along which the switching function 
vanishes identically, see Definition 3.5 below) is always contained in the set A5!1(0). 


Then, the synthesis is built recursively w.r.t. the number of switchings of extremal trajectories, 
canceling at each step the non optimal trajectories (see [12, Chapter 1]). 
In the optimal syntheses some special curves appears, namely: 


e Switching curves, i.e., curves made of points where the control switches from +1 to —1 
Or vice versa. 


e Singular curves, i.e., curves along which the corresponding switching function (see Defi- 
nition 3.5 below) is identically zero. For these curves the corresponding control can take 
values in the interior of [0,1]. 


e Cut loci, namely curves made of cut points. Given an admissible curve 7 of the control 
system, defined on [0,7], we say that teut € (0,7) (resp. 7 (teut)) is à cut time (resp. 
point) if Jo] is time-optimal and y|pr+e] is not for every sufficiently small €. 


Remark 3.4. Notice that, although the functions À 4 and AB depend on the coordinate chart, 
the sets A3'(0), AE (0) and the function fs do not, i.e., they are intrinsic objects related with 
the control system (3.3). 


We are now interested in determining the extremals. À key role is played by the following: 
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Definition 3.5 (switching function). Let (q(-),p(-)) be an extremal pair. The corresponding 
switching function is defined as @(t) = (p(t), G(q(t))). 


Notice that @(-) is continuously differentiable since @(t) = (p(t),[F,G](q(t))) which is con- 
timuous. 


Definition 3.6 (bang, singular). Let q(-) be an extremal trajectory defined on the time interval 
[a,b], and let u(-) : [a,b] — [-1,1] be the corresponding control. We say that u(-) is a bang 
control if u(t) = +1 a.e. in [ab] or u(t) = —1 a.e. in [a,b]. We say that u(-) is singular if 
the corresponding switching function @(t) = 0 in [a,b]. A finite concatenation of bang controls 
is called à bang-bang control. À switching time of u(-) is a time 7 € {a, b] such that, for every 
€ > 0, u(-) is not bang or singular on (r—€,r+€)f{[a,b]. An extremal trajectory of the control 
system (3.3) is called à bang extremal, singular extremal, bang-bang extremal respectively, if 
it corresponds to a bang control, singular control, bang-bang control respectively. If 7 is a 
switching time, the corresponding point q(r) on the trajectory q(-) is called a switching point. 


The switching function is important because it determines where the controls may switch. 
In fact, using the PMP, one easily gets: 


Proposition 3.7 ([12, Lemma 5 page 40]). À necessary condition for a time 7 to be a switching 
time is that O(T) = 0. Therefore, on any interval where has no zeros (respectively finitely 
many zeros), the corresponding control is bang (respectively bang-bang). In particular, 9 > 0 
(resp. d < 0) on [a, b] implies u = 1 (resp. u = —1) a.e. on [ab]. On the other hand, if d has 
a zero at T and d(T) is different from zero, then T is an isolated switching. 


8.2.2 More on singular extremals and predicting switchings for 2D-systems 


In this section we compute the control corresponding to singular extremals and we would like 
to predict which kind of switchings can happen, using properties of the vector fields F and G. 
The following lemmas illustrate the role of the set AE (0) in relation with singular extremals. 
Proofs can be found in [11, 12, 36|. 

We denote by Supp (q(-)) the support of the trajectory q(-). 


Lemma 3.8 ([12, Lemma 10 page 46[). Let q(-) be an extremal trajectory which is singular on 
a time interval [a, b] included in the domain of q(-). Then, q(-)law corresponds to the so-called 
singular control w(q(t)), where 

d AB F q 

CT Ge 

daAp (G(a)) 
with AB defined in equation (3.6) and d,Ag being the differential map of Ag at point q. 
Moreover, on Supp (q(-)), w(q) is always well-defined and its absolute value is smaller or equal 
to one. Finally Supp (q(-)lab) € NEO): 


(3.8) 


The following lemma describes what happens when A1 and A3 are different from zero. 


Lemma 3.9 ([12, Theorem 11 page 44]). Let Q C M be an open connected set such that 
QAN(A;!(0)UA;1(0)) = 9. Then all connected components of Supp(q(-))NQ, where q(-) is 
an extremal trajectory of (8.8), are bang-bang with at most one switching. Moreover, if fsla > 0, 
then q(-)la is associated to à constant control equal to +1 or —1 or has a switching from —1 
to +1. 1f fsla < 0, then q(-)la is associated to a constant control equal to +1 or —1 or has a 
switching from +1 to —1. 
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Remark 3.10. For the problem (R), under generic conditions on the vector fields F and G, 
one can make the complete classification of synthesis singularities, stable synthesis, singularities 
of the minimum time wave fronts. We refer to [12] for the general theory. For some extensions 
to higher dimension, see [2, 41]. 


3.3 Construction of the time-optimal stabilizing synthesis for the reduced 
system 


Let us apply the theory of the previous section to problem (Q). 


3.3.1 Preliminary observations 


Define 


r-(s). a(#), v=r+e= (it), x=r-@- (rt). 
0 — —T Fa 


The integral curves of the vector fields Y and X are shown below (Fig. 4). 


é LL 2 LR NN 


\ 7 
T CL | 
RS 


3 2 -1 0 1 2 3 


Figure 4: Integral curves of X and Y 


We compute the quantities, 
1 ÿ à 
AA (&,9) = det(F(&,ÿ), G(&,ÿ)) = der (* :) = $; 


Le Là y 0 2 
Ao(,D) = det(G(&. D) [F,GI(E D) = det | 3, 9) = à 
e Ag(&,ÿ) ÿ 
Ë,9) = === = 2. 
Js(E,9) AD 
Using Lemmas 3.8 and 3.9 it follows that: 


e on the set {& > 0} f{ÿ > 0} only the switching from +1 to —1 is allowed. And cyclically 
on the other quadrants. 
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e If there are singular trajectories then they should lie on the set {ÿ = 0}. Notice that a 
trajectory whose support is {ÿ = 0} ? is admissible (it corresponds to the zero control). 
See Fig. 5. 


fs > 0 fs <0 


fs <0 fs >0 


Figure 5: Candidate singular 


We construct the time-optimal synthesis in the following way: 


e STEP 0. We consider the trajectory 7 starting from qo with control —1. All extremal 
trajectories should start in this way since in a small neighborhood V of qo we have that 
fs < 0 and hence neither singular trajectories nor trajectories having chattering, i.e. for 
instance an infinite number of switchings in finite time, could be extremal in V. Moreover, 
a trajectory cannot start from qo with control +1 since the vector field Y (go) = 0. 


e STEP 1. We compute the last time at which the trajectory + is extremal (or has lost 
optimality because it self-intersects) and we study which kind of extremal trajectories 
can bifurcate from 7. 


e STEP 2. For each bang or singular trajectory bifurcating from y, we study the last time 
at which it is extremal. If there are intersections among trajectories, we cancel those 
parts that are not optimal (among trajectories already computed up to this step). 


e STEP 3. For each trajectory 7 defined on [0,7] computed at the previous step do the 
following. If at 7 this trajectory loses optimality (among trajectories of the previous 
step), do nothing otherwise prolong + with the next bang or singular trajectory up to 
the last time at which it is extremal. If there are intersections among trajectories, cancel 
those parts that are not optimal (among trajectories already computed up to this step). 
Notice that in general there could be several prolongations of + that are extremals (due 
to the presence of singular trajectories, see [12, $2.6.4 page 62]). Hence, in general, new 
extremal trajectories are generated at this step. 


Then, the synthesis is built recursively, repeating STEP 3 up to when no new trajectories are 
generated. Notice that a trajectory generated at STEP n has n concatenations of bang or 
singular trajectories. In our case, only one application of STEP 3 is necessary. 


?According to [12], the set {ÿ = 0} is a turnpike (i.e., a trajectory that is locally time-optimal) 
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8.3.2 Adjoint equation 


To compute the last time at which bang and singular trajectories are extremal, we need the 
adjoint equation. From equation (3.4) we have that H(&,ÿ,£,0,u) = —£6 + u(£ÿ — C&) (here 
q = (&,ÿ) and p = (£,6)). Hence, , 

{ É = u(t)c 


: 3.9 
= -u(té. + 


3.3.3 STEP 1: Analysis along the curve starting from 4 and corresponding to 
control —-1 


The trajectory y starting from qo and corresponding to control —1 has coordinates: 


This curve intersects the &-axis at time {1 — 7/3. The solution of the adjoint equation with 
u = —1 and the normalization |p(0)| = 1 is, 


E(t, a) = cost + à) 
G(t,a) = sin(t + a), 


where & is defined by cos(a) — £(0), sin(a) = €(0). Notice that the condition H + À = 
E + u(éÿ — CE) + À = 0 with À < 0, written at the initial point implies that, 


— cos(a) + ucos(a) = cos(a)(u — 1) > 0. 


Since we start with control —1 we have that a € [r/2,3/2r|. 
Then, the switching function is given by 


bat) = (p(#), G(17 (D) = 2cos(a) — cos(t + a). 
By studying this function it follows that: 


e for a = x/2, 4, starts with value zero and then takes positive values: it cannot correspond 
to à trajectory starting with control —1. 


e Forae (x/2,2r/3), pa starts with negative values, and has its first zero (of order one) 
(close to the origin for values of a close to x/2). For every Tr € (0,7/3) there exists a 
unique @ € (x/2,7/3) for which the switching function has a zero at 7. 


e For à — 27/3, da starts with negative values and has à double zero at t = æ/3 (the 
point gs at which the trajectory meets the turnpike). Notice that there exists an ex- 
tremal trajectory entering the turnpike, since there is a trajectory switching there (see 
the singularity (Y, S) in [12, page 62]). 


e For a € (2r/3,47/3), Pa is always negative. This means that the trajectory with control 
—1 is always extremal (a priori, y can make many turns !). 


e For a = 47/3, da starts with negative values and has à double zero at t = 57/3 (again 
at à point gs where the trajectory meets the turnpike). Again notice that there exists an 
extremal trajectory entering the turnpike, since there is a trajectory switching there. 
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e Forae (47/3,37n/2), ÿ starts with negative values, and then has its first zero (of order 
one) (close to t = 57/3 for small values of a close 4/3). For every 7 € (x,57/3) there 
exists à unique @ € (4x/3,37/2) for which the switching function has à zero at 7. 


e For à — 37/2, ®4 starts with negative values, and has its first zero (of order one) at 
É= TT. 


This study being done, we get the expression of the first switching time 71 solving the equation 
Palri) = 0 (and taking into account that 71 > 0). It yields 


na) =-a+arccos(2cosa), a€ (5 æ U Es = : (3.10) 


Notice that y is extremal for every t > 0. However, y cannot be optimal after t = 27 since 
it is periodic. 
The initial conditions of the covector are shown below (Fig. 6). 


édité tE (0,7), pa has one zero 


— , _1p2(0) 
ET x/2 
Ba € Ô 4 
IT p1(0) 
| En 
\ 17 
 — 2 
nes ed oi 


Kb every t € (Tr, 2x), da has one zero 


Figure 6: Initial conditions of covector 


3.3.4 STEP 2 
At this step we study separately the following four families of trajectories bifurcating from 7". 
e Family 1: trajectories bifurcating from +7 [0,x/3). 
e Family 2: the singular trajectory bifurcating from gs = y-(x/3) = (—-V3,0). 
e Family 3: trajectories bifurcating from +7 [r,5/3r] 
e Family 4: the singular trajectory bifurcating from gs = +" (5x/3) = (13,0) 


The extremal trajectories bifurcating from + are drawn in the following picture (Fig. 7). 


19 


71 


Figure 7: Extremal trajectories bifurcating from y 


Family 1. Consider a curve corresponding to control +1 starting from the point (71), with 
n = (a), à € [r/2,2x/3] (which corresponds to 1 € [0,w/3]). Its coordinates are: 


x(t, 7) = 2sin(t - nn) — 2sin(t) 
y(t, T7) = 1—2cos(t) + 2cos(t — T1) 


The corresponding covector which satisfies the switching condition at the beginning 
p(0, n)G(x(0, Ti); y(0, T1)) _ 0, 
has coordinates 


E(t, 71) = cos(t — 71 — à) 
GE, 71) = —sin(t — 7 — à), 


The switching function is given by (here we have taken account that d,(r1) = 0) 


Palt + Ti) = cos(t — T1 — à) — 2 cos(a). 


We get the expression of the second switching time 72 solving the equation Pal +7) = 0 
(and taking into account that r > 0). It yields 


(3.11) 


2 
m(a) =2arccos(2cosa), € (£ e| . 
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The switching curve is shown in the following picture (Fig. 8). 
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switching curve 


Figure 8: Switching curve family 1 


Family 2: the singular trajectory. Since there is a trajectory switching from y in qs 
it follows that there exists an extremal trajectory entering the turnpike. At the entrance 
point of the turnpike the covector satisfes (ps, G(qs)) = 0. Since G = (y, —-x) it follows that 
(ps, (0, V3)) = 0. Hence we can assume that 


Ps — A0): 


The minus sign has been chosen in such a way that (ps, F + uG) > 0 (according to Theorem 
3.1, H+À= (ps, F +uG) + À = 0, with À < 0). On the singular trajectory the control can 
be only zero. From equation (3.9) it follows that p is constant. Hence the singular trajectory 
is extremal for all positive times (Fig. 9). 


Figure 9: Singular arc 


Family 3. This family is treated exactly as Family 1. The only difference is that now à € 
Ar/3,37/2] (which corresponds to 7 € [r,5x/3]). The switching curve is shown on the 
following picture (Fig. 10). 
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Figure 10: Switching curve family 3 


Family 4. This family is treated exactly as Family 2. Since there is a trajectory switching 
from y at gs it follows that there exists an extremal trajectory entering the turnpike. This 
trajectory is then extremal for every positive time (Fig. 11). 


ÿ 
q 


Figure 11 


Cutting non optimal trajectories. À direct computation yields the following. 


e The trajectory y is not optimal after passing through the point gs. 


e Let us Call Jcrazy the trajectory of Family 3 corresponding to r1 = 7 (i.e., & = 37/2). 
This trajectory is not optimal after crossing the singular trajectory of Family 2. 


e The trajectories of Family 3 for n € (x,57/3) (i.e., a € (41/3,37/2)) 
are not optimal after crossing the 77 trajectory. 


e Let us call ÿaoubie the trajectory of Family 3 corresponding to r1 = 57/3 (1.e., « = 4x/3). 
This trajectory is not optimal after passing through the point gs. Notice that the support 
of this trajectory coincides with part of the support of à trajectory of Family 1 (the one 
corresponding to 71 = 7/3 (i.e., à = 27/3). 
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e The singular trajectory of Family 4 is not optimal. 


It follows that only the trajectories shown in the following picture are candidates to be optimal. 
These are the trajectories generated at STEP 2 (Fig. 12). 


Y 


Figure 12: Trajectories generated at STEP 2 


3.3.5 STEP 3 


The trajectories generated at this step are shown on the following figure where they are divided 
in four families: family 2.1, 2.2, 2.3 and 2.4 (Fig. 13). 


y 


family 2.1 


Lit 


Figure 13: Trajectories generated at STEP 3 


The point qtop that divides the families 2.3 and 2.4 is the point at which the switching curve 
generated at the previous step become tangent to an integral curve of the vector field X. 


The coordinates of the point Gtop are ( — 1/ Ey3 — 51, à (5V3 — 9) ). 


Remark 3.11. Notice the following important fact. The trajectory of Family 1 of STEP 1 
corresponding to a = arccos (—5/4V2) (i.e., r1 = 2arcsin(1/4) and > = 37/2 + arcsin(1/4)) 
has a switching at the origin at time { = 37/2 +3 arcsin(1/4). However, this trajectory cannot 
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be connected to a singular trajectory. This is due to the fact that the vector field G is zero at 


this point. Hence, the switching function vanishes but the covector has not necessary the right 
direction to make the trajectory enter the turnpike extremal. This phenomenon is similar to 
the one described in [10, Fig. 4.6]. A direct computation shows that the value of the covector of 
this trajectory at the origin is (— V5/8, 3/3) while à covector of the form (—1,0) is necessary 
to enter the turnpike. 


We have the following 


The family 2.4 is clearly not optimal since it is generated by a switching curve that reflects 
the trajectories (see the analysis of C curves in [12, $4.1.3/). 


The families 2.1 and 2.2 do not stop to be extremal before meeting the set K1 = {(x,0)|x > 
V3}. This set is a cut locus. The curve crazy is not optimal. It is quicker to take 
trajectories from the family 2.2. 


The family 2.3 loses optimality when it meets Family 3 of the previous step. It forms then 
a cut locus 2 which is à one-dimensional manifold starting from gp and arriving to a 
point Gbottom On the y trajectory. The coordinates of the point {bottom are (—V3, —2). 


This cut locus has the following features: 


— At Gtop it is C l_connected to the switching curve from which the family 2.3 is gener- 


ated. This tangency is a consequence of the general theory developed in [12]. It is 
called à (C, K)2 singularity (see [12, Fig. 2.21]). 


It is not smooth at one point that we denote by Qmiaae: More precisely, it is not 
smooth at the intersection with the curve Jaoubie that separates the families 1 and 3 
of the previous STEP. The presence of this nonsmoothness point is a consequence of 
the fact that the minimum-time front is not smooth along the curve oublie. Indeed 
this trajectory separates two different families. 


The coordinates of the point Qidatle are (-Æ, -à). 


At the point dbottom it is C'!-connected to the curve 7. This fact is the consequence 
of the tangency of Family 3 with y at the point y” (x). 


Since all the trajectories generated at this step are either defined on [0, +) or on [0,7] with 
cut time T,, it follows that no new trajectories are generated at subsequent steps. Hence, the 
time-optimal synthesis is constructed (Fig. 14). 


Remark 3.12. Notice that the minimum time function is not continuous on +7 ([0,x)). 
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Figure 14: The optimal synthesis. All optimal trajectories start at go with control —1. The 
dashed black curve is the switching curve, the purple curve is the singular trajectory and the 
green curves are cut loci. Notice that the minimum-time function is not continuous along the 
blue dashed curve. 


3.4 Correspondance between the solutions to problem (P) and problem (Q) 


The solutions to problem (P) can be deduced from the solutions to problem (Q). In this 
section we display pairs of figures showing a solution of problem (P) and its corresponding 
lifted solution to problem (Q). Pictures are shown for each type of solution: bang-bang, 
bang-bang-bang, and bang-singular-bang trajectories. 

Let us make a few remark on the time-optimal synthesis of problem (P). 

As the reader can easily check using a similar argument as it was done in [44], the application 
of PMP to problem (P) leads to the following lemma describing some of the feature of the 
synthesis. 


Lemma 3.13. Optimal paths solving Problem (P) are such that rectilinear segments and points 
of inflection belong to a same line passing through the center of the final manifold C. 
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Figure 15: À bang-bang-bang optimal trajectory. À bang-bang-bang solution to problem (Q) 
(left) and the corresponding optimal solution to problem (P) (right). 
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Figure 16: À bang-bang optimal trajectory. À bang-bang solution to problem (Q) (left) and 
the corresponding optimal solution to problem (P) (right). 


4 Tracking a moving target 


The steering methods presented in the previous sections of this paper were concerned with a 
fixed target. In this section, the results of the previous sections are brought together for the 
development of a control strategy for tracking a moving target. 

This tracking problem has been tackled in different ways. For instance, while some authors 
use predefined pieces of trajectories called patterns that depend on the difference of velocity 
between the target and the UAV [30, 34] some others prefer to use some simple algorithm 
based on tangent vector or lateral guidance laws [39, 47] in order to compute in real time the 
desired path for the UAV. Some authors consider that the tracking problem can be modeled 
by a predator-prey model which gives waypoints to be interpolated by the UAV [4]. 

Also, some authors use the same kind of methods as the ones we have described in the 
previous sections. For instance, in [33] a method using “Lyapunov vector fields” is presented to 
track à moving target. 

Time-optimal problem to manage convoy protection has been considered in [21] where the 
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Figure 17: Another bang-bang optimal trajectory. À bang-bang solution to problem (Q) (left) 
and the corresponding optimal solution to problem (P) (right). 
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Figure 18: À bang-singular-bang optimal trajectory. À bang-singular-bang solution to problem 
(Q) (left) and the corresponding optimal solution to problem (P) (right). 


authors consider that the minimum turning radius of UAV is larger than the radius of the 
protection circle around the target, ï.e., the UAV cannot follow this circular trajectory. 


4.1 Statement of the problem 


Our goal is to make the UAV reach and track a time-parameterized circular final manifold C; 
lying in an horizontal plane above the moving target. Both cases when the future trajectory of 
the target is known or not are considered. 

Let 2 — (2,22,z3) represent the target coordinate in R%. Define C; to be the time- 
parameterized counterclockwise oriented circle with minimal turning radius r = 1/umax Cen- 
tered at (z1(6), z2(t)) in a horizontal plane of a given altitude above the target. 

The simple algorithm presented below depicts an open-loop controller which can be used in 
a high-level controller that is itself used to find waypoints to be sent to the low-level controller 
performing the closed-loop control of the UAV. 

Here for every position of a HALE or MALE UAV, modeled by (2.1) or (2.5) respectively, 
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we would like to find a trajectory steering the considered UAV to the minimum curvature hor- 
izontal circle centered at the centroid of a moving target, at a given altitude w.r.t. the target. 


It is necessary to make some assumptions on the target behavior. First, we assume that 
the target speed vtarget is smaller than the UAV speed vuav, since it is clear that without this 
assumption the UAV cannot follow the target. 


Summing up we consider the following problem 


Problem T Given two initial conditions 2(0) = z0 and q(0) = go in the corresponding state 
spaces of the target and the considered UAV, find a pair trajectory-control joining go to Cx which 
is admissible for the considered control system. 


4.2 Description of the algorithm 


We assume that the past behaviors of both the UAV and the target are known. In practice, 
the target trajectory is obtained via an external module (as for example a vision based module 
see [23]). 


The following rough algorithm describes a high-level open-loop target tracking controller. 
We construct the path in the following way 


e STEP 0. Acquisition of initial data z9 € R° and qo € R° x SO(3), the position of the 
target and UAV respectively. 


Set t = 0. 
e STEP 1. Compute t" = min{seR, | q(t+s) e C)}. 
e STEP 2. Predict the value of z(t + t*). 


e STEP 3. Compute the pair trajectory-control (q(-),u(-)) to reach the C;:# and apply 
the control to the UAV. 


e STEP 4. Set t = t + At and go to STEP 1, where At is some time step, that may either 
be constant or updated in terms of external informations. 


Remark 4.1. There is of course no convergence result for this strategy. Here, we just give an 
idea on how to apply our previous steering methods in the case of a moving target. Of course 
the positioning error occurring at STEP 3 of the previous algorithm should be studied. 


Remark 4.2 (more precisions on the different steps of the algorithm). 

-At STEP 1 we compute the final time t* to reach the C;. Depending on the method, the 
value of this final time will be obtained exactly or not. The time-optimal method yields the 
exact time to reach the final manifold C;. In this case, t* is solution to problem (P) (or (Q), see 
Section 3). When the Lyapunov method is used, C; is reached asymptotically and the time t* is 
the first positive time at which the UAV enters a neiborhood C? = B(z(t),r+e){ B(z(t),r —-€) 
of C;. The parameters € can then be tuned in order to improve on the performance of the path 
following algorithm. 


-At STEP 2 we make a prediction of the future position (in space) z(t + t*) of the target. 


Of course this position may be obtained in different ways. The method we have used consists 
of using an observer in order to reconstruct the target velocities. In this work we implemented 
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the so-called extended Kalman filter (EKF for short) (see e.g. [6, 42, 38]) which is well-known 
for its robustness w.r.t. the noise. 
For this purpose, let us assume that the target dynamics is also modeled by an extended 

Dubins system (see [18]) 

1 = U1 cos Ô 

22 = u sin 
“ (4.1) 
23 = U2 
Ô — U3; 


where the controls u1, u2 and u3 are the target horizontal linear velocity, altitude and yaw 
angular velocity respectively. The available measurements are the cartesian target coordinate 
21,22, Z3, thus, according to (4.1), ui, u2,u3 are clearly observable (we used the local model 
ü = 0, i = 1,2,3). Note that system (4.1) reduces to the Dubins system when the target 
moves in à horizontal plane. 


4.3 Numerical results 


We now present numerical simulations illustrating the tracking algorithm presented in the 
previous section. 

In each case we consider a UAV flying at constant unitary speed and modeled by System 
(2.1) or System (2.5) depending on the type of UAV we consider (HALE or MALE). We 
numerically integrate the system during 500 units of time. 

For each of the three methods we show simulations for both cases where the target path is 
known or not. If this information is unknown we use an extended Kalman filter to predict the 
future target position. Note that the target speed is not assumed to be constant. 

On each figure, the dashed circle represents the trajectory to be followed. The results are 
drawn in the inertial frame (on the left) and in an orthonormal moving frame attached to the 
target (on the right) in order to observe the positioning error. 

Figures 19 and 20 illustrate the time-optimal tracking method for the same UAV starting 
from the same initial position. The target that we want to track starts at (21,22) — (0,0); its 
behavior is known on Figure 19 and unknown on Figure 20. 

The second couple of figures, Figures 21 and 22, illustrate the Lyapunov tracking method 
for the same UAV starting from the same initial position. The target that we want to track 
starts at (&o, ÿo) — (0,0); its behavior is known on Figure 21 and unknown on Figure 22. 

Figure 23 and Figures 24, illustrate the 3D-Lyapunov tracking method for the same MALE 
UAV starts horizontally from the point (10, —20, 50) with a yaw direction of 0 radian and that 
has the curvature constraints U1 max = U2max = 0.1, U3max = 0.5, a3 = 0.1. The values of the 
parameters € and © (see Section 2.2) are 1.2 and 1 respectively. The target that we want to 
track starts at (0,0,0) and is supposed to move in an horizontal plane: its behavior is known 
on Figure 23 and unknown on Figure 24. In each simulation, C; lies in the horizontal plane 
located at an altitude of 10 units over the target. 

Finally, Figure 25 shows a tracking trajectory obtained for a MALE UAV which follows a 
target which is not moving in a plane. Here the drone must be positioned at 10 units above 
the target. 
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Figure 19: Time-optimal based tracking method in 2D. UAV'Ss trajectory is continuous, target’s 
trajectory is dashed and the C; is dashed. The behavior of the target is known. (x9, Yo, 00) = 
(15,10,0), (21,22) = (0,0), and wmax = —@ = 0.5. (a) Trajectories drawn in the (x,y)- 
coordinates plane. (b) Trajectories drawn in a (orthonormal) frame attached to the center of 
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Figure 20: Time-optimal based tracking method in 2D. The only difference with Figure 19 is 
that, on this figure, the behavior of the target is unknown. 
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CHAPITRE 3. RÉSUMÉS DES CONTRIBUTIONS PERSONNELLES 


3.3. Généralisation des méthodes pour atteindre n’importe quel pattern 


3.3 Généralisation des méthodes pour atteindre n’importe quel pat- 
tern 


Dans le paragraphe précédent, nous avons proposé des méthodes de planification permettant de 
rejoindre un cercle. 

Dans le but de planifier une trajectoire, pour tous types de missions, il faut que le pattern visé 
puisse être différent d’un cercle (e.g., un hippodrome ou un huit), cf. paragraphe 1.1. 


En s'inspirant de [81], pour rejoindre n'importe quel pattern, nous utilisons une des méthodes 
décrites précédemment en visant un cercle inscrit au pattern considéré. 

Nous utilisons deux cercles inscrits pour définir les patterns en forme d’hippodrome et de huit. 
La position de ces cercles définit l’aspect du pattern (e.g., longueur, orientation). 

Un exemple d'application est représenté en Figure 3.8. Pour calculer ces trajectoires, nous avons 
utilisé la méthode de planification temps-minimal. 


FIGURE 3.8 — Une trajectoire optimale en temps reliant (a) un hippodrome, (b) un huit. 


3.4 Comment l’appareil peut suivre une cible dans un milieu en- 
combré ? 
Le but est d'appliquer les résultats exposés précédemment pour effectuer du suivi de convoi en 
milieu encombré. Cette planification est possible en s’inspirant de : 
1. La planification point-point en milieu encombré (cf. paragraphe 2.5.3), 
2. La méthode de suivi de convoi en milieu non contraint (cf. paragraphe 3.2). 


Les notations sont celles du paragraphe 3.2. La trajectoire est construite de la manière suivante : 
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CHAPITRE 3. RÉSUMÉS DES CONTRIBUTIONS PERSONNELLES 


3.4. Comment l’appareil peut suivre une cible dans un milieu encombré ? 


Étape O0 Acquisition des données initiales cible) € R° et go € R° x SO(3), les positions de la 
cible et du vecteur respectivement. 
Poser t := 0. 


Étape 1 Calcule de t* = minfr ER, | q(t+r) e &)}. 
Étape 2 Prédiction de la valeur de qibie(t + t*). 


Étape 3 Calcule du couple (q(-),u(-)) permettant de rejoindre le pattern €,# en utilisant la 

méthode suivante : 

Étape 3.0 Initialisation du point courant (position du vecteur). 

Étape 3.1 Recherche du fil d'Ariane reliant le point courant à dcible(t + t*). 

Étape 3.2 Recherche d'une trajectoire reliant le point courant et un pattern centré autour 
d’un point du fil d'Ariane et le plus éloigné possible. 

Étape 3.3 Mise à jour du point courant et retour à l'étape 3.1, tant que celui-ci n’est pas 
dans un voisinage du point d’arrivée. 

Étape 3.4 Appliquer le contrôle au vecteur. 


Étape 4 Poser t := {+ At et recommencer à partir de l'étape 1. At est un pas de temps fixe ou 
actualisé en fonction d’un module externe. 


Remarque 3.6. Le calcul de t* à l’étape 1 peut être obtenu en appliquant une procédure similaire à 
celle de l'étape 3. 

La méthode de planification point-pattern de l’étape 3.2 est choisie parmi celles présentées au 
paragraphe 3.2. 


Des exemples d'applications sont donnés en Figures 3.9 et 3.10. La trajectoire du convoi suivi est 
en pointillés. Une méthode de planification point-pattern temps-minimal a été utilisée pour obtenir 
la trajectoire de la première figure (Fig. 3.9). Nous avons utilisé la méthode basée sur le principe de 
LaSalle pour calculer la trajectoire de la Figure 3.10. 
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3.4. Comment l’appareil peut suivre une cible dans un milieu encombré ? 


FIGURE 3.9 — Planification de trajectoires en milieu contraint utilisant la méthode temps-minimal. 
La trajectoire du convoi est en pointillés. 
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FIGURE 3.10 — Planification de trajectoires en milieu contraint utilisant la méthode basée sur le 
principe de LaSalle 2D. La trajectoire du convoi est en pointillés. 
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3.4. Comment l’appareil peut suivre une cible dans un milieu encombré ? 
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Chapitre 4 


Conclusion 


Le travail effectué dans cette première partie concerne la planification de trajectoires pour des 
drones de type HALE ou MALE, qui peuvent être modélisés par des systèmes différentiels non- 
linéaires simples. 


Le chapitre 2 nous a permis d'exposer les méthodes usuelles de résolution de ce type de problème. 
Les premiers résultats que nous avons obtenus pour la planification vers un point, ont été étendus 
r) 
pour prendre en compte des contraintes environnementales. 


Le chapitre 3 a apporté des réponses au problème de couplage vecteur/capteur. Nous avons 
proposé une méthode d’asservissement de la trajectoire du couple drone/caméra donnant plus ou 
moins d'importance à l’un ou à l’autre par le biais de pondérations. 

Nous avons remarqué, dans le cas d’une caméra non contrainte, que ce problème de couplage se 
réduit au problème de planification du vecteur. Dans ce cas, une méthode indépendante de contrôle 
de la caméra a été proposée. 

Le problème de planification a été traité pour une mission impliquant l’utilisation de patterns. 
Nous avons exposé deux méthodes de planification point-pattern. Ces deux méthodes, basées res- 
pectivement sur le principe de LaSalle et sur le PMP, ont été testées numériquement. 

Pour finir, nous avons proposé et appliqué un algorithme de suivi de cible, basé sur les méthodes 
de planification vers un pattern. 


Dans toute cette première partie, des algorithmes de planification point-point et point-pattern 
ont été étudiés. Parmi toutes ces méthodes de calcul de trajectoires, certaines sont obtenues par la 


minimisation d’une fonction de coût. 


Dans la partie suivante, nous proposons une méthode de construction d’une fonction de coût 
permettant de calculer des trajectoires les plus proches possibles de celles effectuées par les pilotes. 
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Deuxième partie 


Une méthode anthropomorphique : 
l’utilisation du contrôle optimal inverse 
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Chapitre 5 


Introduction 


Dans la partie précédente, nous avons étudié des méthodes de planification point-point et point- 
pattern. Parmi toutes ces méthodes de calcul de trajectoires, certaines sont obtenues par la minimi- 
sation d’une fonction de coût. 

Le problème que nous allons étudier, dans cette partie, est le suivant : comment proposer des 
trajectoires voisines de celles qu’un pilote expérimenté serait amené à décrire s’il était directement 
aux commandes de l’appareil ? 

En d’autres termes, à partir du comportement d’un aéronef, nous souhaitons extraire un coût 
afin de l’utiliser au sein des modules de planification précédemment présentés. 


Cette approche, connue sous le nom de contrôle optimal inverse, a été appliquée à des problèmes 
de biomécanique humaine (quelques références sont données page 101). Cette méthode d'analyse de 
données très novatrice est basée sur un paradigme aujourd’hui largement accepté en neurophysiologie 
selon lequel, parmi tous les mouvements possibles ceux accomplis satisfont un critère d’optimalité 
[134]. 


L'utilisation de cette analyse est importante puisqu'elle nous permet d’appliquer les règles de 
l'Organisation de l'Aviation Civile Internationale, et notamment l’article 2.13 de la circulaire 328 de 


l'OACI [74]. 


Article 2.13 de la circulaire 328 de l'OACI 


À key factor in safely integrating UAS in non-segregated airspace will be their ability to act and 
respond as manned aircraft do. |... 


International Civil Aviation Organization 


Dans le chapitre 6, nous exposons quelques travaux déjà effectués dans le domaine. Le chapitre 
7 présente nos propres résultats qui constituent en fait un article accepté dans la revue ESAIM : 
Control, Optimisation and Calculus of Variations (à paraître en 2013). 
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Chapitre 6 


Méthodes bio-inspirées et problème de 
contrôle optimal inverse 


Sommaire 
6.1 Les méthodes bio-inspirées ........................... 100 
6.2 Le problème de contrôle optimal inverse .................. 100 
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CHAPITRE 6. MÉTHODES BIO-INSPIRÉES ET PROBLÈME DE CONTRÔLE OPTIMAL 
INVERSE 6.1. Les méthodes bio-inspirées 


6.1 Les méthodes bio-inspirées 


La bio-inspiration, ou biomimétisme, est une philosophie de développement basée sur l’inspiration 
par la nature. 

De nombreux travaux ont proposé de modéliser et créer des robots en mimant des insectes ou 
des animaux. C’est le cas par exemple dans [112] où les auteurs s’inspirent du comportement des 
mouches pour effectuer de la planification sommaire de trajectoires, afin de longer un mur. 

D’autres encore s’inspirent des criquets pour créer des capteurs de positionnement appliqués à 
des véhicules aériens [114]. 

Certains tentent d’imiter le comportement des jongleurs pour effectuer de la manipulation d’ob- 
jets avec des bras robotiques [7]. Cette même équipe à aussi travaillé sur un robot mimant une 
salamandre [65]. 


Tous ces travaux ont été possibles grâce à des modélisations et des études de comportements na- 
turels. Par exemple, le corps humain a beaucoup été étudié [67, 118, 9, 11]. De plus, grâce à des études 
comportementales, certains principes naturels ont été découverts, e.g., le principe d’inactivation [20]. 


Afin de traiter ce problème de biomimétisme, certains auteurs considèrent que la nature peut 
être analysée via un problème de contrôle optimal inverse [30, 87, 17]. Ce problème est explicité dans 
le paragraphe suivant. 


6.2 Le problème de contrôle optimal inverse 


Le problème de contrôle optimal inverse est similaire à de la cryptanalyse. Rappelons d’abord 
que la cryptanalyse a pour but de retrouver la clef permettant de décoder un texte chiffré, écrit 
par une tierce personne. Le principe de résolution du problème de contrôle optimal est le même : 
il permet de trouver le critère d’optimalité (la clef) d’un problème de contrôle optimal, à partir 
d'observations (le texte à décrypter). 


Comme pour effectuer une cryptanalyse, nous avons besoin de poser des hypothèses plausibles 
sur le phénomène étudié. Ces hypothèses concernent en général la classe de systèmes différentiels 
choisie pour modéliser le phénomène, ainsi que la classe de régularité de ces systèmes et du critère 
recherché. 

Il faut ensuite retrouver, à partir d'observations, le critère d’optimalité minimisé. Le problème 
se formule de la manière suivante 


Problème de contrôle optimal inverse 


Étant donnés un système de contrôle 4 — f(q,u) et un ensemble T de trajectoires admissibles 
observées, trouver un critère C(q(-),u(-)) tel que toute trajectoire 7 € T° soit solution de 


inf{C(q(},u()) | 4= fau), g(0)=7(0), a(T) = 7(7)} 
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CHAPITRE 6. MÉTHODES BIO-INSPIRÉES ET PROBLÈME DE CONTRÔLE OPTIMAL 
INVERSE 6.2. Le problème de contrôle optimal inverse 


En acceptant le paradigme de neurophysiologie selon lequel, parmi tous les mouvements possibles 
ceux accomplis satisfont un critère d’optimalité [134], les actions humaines sont analysables via un 
problème de contrôle optimal inverse. 
Voici une liste, non exhaustive, des phénomènes étudiés : 
e Mouvement de pointage du bras (B. Berret, J.-P. Gauthier, F. Jean, E. Todorov, F. Jean) 
[20, 99, 19, 99] 

e Mouvement de saccade des yeux (F. Jean) 

e Locomotion humaine (Y. Chitour, F. Chittaro, F. Jean, J.-P. Laumond, P. Mason, E. Todorov) 
[10 9: 1L..89, 52, 17, At] 


Une des difficultés, dans ces travaux, est due à l’aspect expérimental de la résolution : la présence 
de bruits de mesure ne facilite pas la reconstruction du coût. Une autre difficulté est la modélisa- 
tion du problème : il n’existe pas toujours de modèle permettant d’expliquer la dynamique d’un 
système biologique. À cause de ces difficultés, la précision du critère reconstruit peut être limitée. 
Néanmoins, des travaux récents, dont celui de F. Jean et al. [11], donnent des résultats de robustesse 
du problème. Ces résultats permettent de savoir si ces incertitudes ont une grande influence sur la 
reconstruction du critère. 


Dans le chapitre suivant, nous modélisons et étudions le problème de contrôle optimal inverse 
pour l’appliquer aux trajectoires de pilotes expérimentés. 


101 


CHAPITRE 6. MÉTHODES BIO-INSPIRÉES ET PROBLÈME DE CONTRÔLE OPTIMAL 
INVERSE 6.2. Le problème de contrôle optimal inverse 


102 


Chapitre 7 


Une méthode d’estimation du coût 
optimisé par des pilotes humains 


L'article présenté dans ce Chapitre a été publié dans le journal ESAIM : Control, Optimisation 
and Calculus of Variations [5]. 
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Abstract. This paper is devoted to the general problem of reconstructing the cost from the observation 
of trajectories, in a problem of optimal control. It is motivated by the following applied problem, 
concerning HALE drones: one would like them to decide by themselves for their trajectories, and 
to behave at least as a good human pilot. This applied question is very similar to the problem of 
determining what is minimized in human locomotion. These starting points are the reasons for the 
particular classes of control systems and of costs under consideration. To summarize, our conclusion is 
that in general, inside these classes, three experiments visiting the same values of the control are needed 
to reconstruct the cost, and two experiments are in general not enough. The method is constructive. 

The proof of these results is mostly based upon the Thom'’s transversality theory. 
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1. INTRODUCTION 


This study holds in the context of the french FUI SHARE project (see [2]), and authors are granted from the 
project. 

From the kinematic point of view, a rough HALE drone (high altitude, long endurance drone) is governed by 
the standard Dubins equations: 


à = cos Ü 
ÿ = sin 0 (1.1) 
Ü=u. 


These equations express that the drone moves on a perfect plane (perfect constant altitude), at perfect constant 
speed 1, moves in the direction of its velocity vector, and is able to turn right and left. 
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These equations were also considered by people studying human locomotion: see the seminal works of Chitour- 
Jean-Mason [8], and Chittaro-Jean-Mason [10]. See also older but basic works of Laumond [3-5], who originally 
introduced this point of view. 

Another typical result related with cost reconstruction in human movements is [6, 7,11]. À part of the 
methodology in [6,7,11] is very similar to what we do here in. 

Note that the rough model (1.1) is more pertinent as a model of the kinematics of HALE drones than as a 
model of human locomotion: actually, for HALE drones, the velocity is almost really constant°. 

Once this assumption of constant velocity is accepted, the natural requirement of invariance under motions 
of the plane leads to these equations. 


In both cases but for different reasons, the assumption is that some integral cost is minimized, to connect in 
free time, the initial and target point. Again, the natural requirement of invariance with respect to motions 
of the plane leads to an integral cost of the form: 


C{u(-)) = 1 L(u(t),üu(t),.…,u(M(t))at (1.2) 


In the case of human locomotion, a very interesting argument of dimension of a certain set shows that k — 1 
in equation (1.2) above (see [13]). We discuss this in more details in Sections 3.3, 5. Shortly, the results here 
provide à clear method to validate this assumption k = 1. 

It turns out, and this is quickly shown in Section 5, that the measurement of two different trajectories is 
enough to reconstruct the cost L(u) (over the set of values of u(t) visited during the two experiments), provided 
that both experiments visit the same value of the control (at different times, maybe). 


These considerations are the starting point of this study, and the reason to consider the following class of 
nonlinear control affine systems, with single control u, and of costs C(u(:-)): 


Le 


ÿ=u 


(1.3) 


where x is n-dimensional, and the cost 


C(u(-)) = Î L(u(t))dt. (14) 


Surprisingly, our main result is the following rough statement, the precise result being Theorem 4.1, 
Section 4.1. 


Theorem 1.1 (rough statement). To reconstruct the cost (1.4), three experiments are in general enough pro- 
vided that they visit the same value of the control, at some (maybe different) instants. Two experiments are in 
general not enough. 


Hence, in fact, the Dubins case is, inside its class, very particular. This is due to the symmetries (invariance 
of the problem under the action of motions of the plane). 

Morcover, in both cases (Dubins case and case (1.3)), the reconstruction of L can be made only up 
to a linear term, which is in fact totally irrelevant: adding a linear term to Z does not change the set 
of optimal trajectories. On the contrary, for completely general affine control systems, this additional degree of 
freedom does not (generically) exist. 


In the paper, we finish by considering (more quickly) the case of à general control affine system, 4 = Æ(q) + 
uF(q) (Sect. 6). It turns out that the main Theorem 1.1 still holds. Even a bit more: now Lu) is completely 
determined by three experiments (not up to some linear term). 


5The argument of possible length reparametrization used in [8, Sect. 2, item (ii), p. 150], seems questionable. 
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The organization of the paper is as follows: 


In Sections 2, 3, we give a clear statement of the problem, together with all prerequisites necessary to the 
mathematical precise understanding of the paper. À key lemma for the proofs (Lem. 3.17) is established. The 
very short Section 3.3 concerns our conclusion about the validation of the assumption that the cost C{u(-)) 
is of the form (1.4). 

- Section 4 is devoted to the results of our study. In Section 4.1 we give a clear statement of our main result 
and in Section 4.2 we give the crucial preliminaries for the proof of the main result (Thm. 4.1). Section 
4.3.2 regroups, besides the final proof, all technical lemmas, and the proof of the fact that in general two 
experiments are not enough. 

- Section 5 deals with Dubins case, and the application to HALE drones. 

Section 6 cares about the case of general control affine systems (very shortly). 


1 


2. STATEMENT OF THE PROBLEM 


In this paper, smooth means C'©. Also, we do not distinguish between row and column vectors, the distinction 
being clear from the context. 


2.1. Systems under consideration 


As we said in the introduction, we deal with single input nonlinear control affine systems of the form (1.3), 
with (x,y) € X x YŸ, the state variable and u € R, the control variable. Here, X (resp. Y) is an n-dimensional 
(resp. one-dimensional) connected, Hausdorff and paracompact smooth manifold. 

Setting q = (x,y) and Q = X x Y, we rewrite the system (1.3) as 


4 = Fo(q) +uFi(q), qe Q, uER, (2.1) 


where F5 = (f,0) and F1 = (0,1) are smooth vector fields on Q. Depending on the context, we will use freely 
the most convenient notation between both. 

Let F denote the set of smooth control systems form (1.3) or (2.1), i.e., elements of F are smooth vector 
fields F5 on Q such that F5 : Q — TX x {0ry}. Equivalently, elements of F are y-parametrized families of 
smooth vector fields f over X, i.e., smooth sections of a vector bundle over Q with fiber T,X. Let JFF denote 
the bundle of k-jets of systems in F. 

The set of admissible control functions u(-) is as usual LX.(R;,R), and the set of admissible trajectories is 
the set of solutions of system (1.3) or (2.1) corresponding to an admissible control function u(-). The domain of 
our control functions will be a certain time interval [f1,t2] depending on u(-). If necessary, it is always possible 
to assume that t1 = (0. 


2.2. Problem under consideration 


2.2.1. Main assumption 
Our main assumption in the paper is: all trajectories occuring in the observations are solutions of an optimal 


control problem of the following form, that we shall denote by (Pr). 


(Pz) Minimize the integral cost (1.4) among all admissible controls u(-) steering system (2.1) from a source point 
do to a target point q1 in free final time T. 


We do not consider arbitrary functions L in the functionals C{u(-)) above: indeed, we choose a class £ 
ensuring that the optimal control problem (P;) has a solution for each L in £. Following [8,10], £ is the set of 
L:R—R meeting the following assumptions: 


(A1) L(-) is strictly positive and smooth; 
(A2) L(:) is strictly convex with L//{u) > 0 (1.e., strictly strongly convex); 


A 
A3) limpy}-+0 L(u)/|u|] = +oo. 
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2.2.2. The different notions of experiments 


Among the set of admissible pairs (q(-),u(-)), we shall distinguish the ones called experiments: 


Definition 2.1 (experiment). An admissible trajectory (q(-),u(-)) is called an experiment if there exists LE L 
such that (q(-),u(-)) solves (Pr). 


Remark 2.2. Since a solution of system (2.1) is uniquely determined by the initial condition qo in Q and the 
admissible control function u(-), we may identify an experiment (q(-),u(-)) with à pair (go, u(-)). For simplicity, 
we shall denote indifferently both by 7. 

We also denote by l; a family of experiments whose corresponding trajectories q(-) are solutions to 
system (2.1). 


Definition 2.3 (compatible experiments). À family 1} of admissible pairs 7 = (q,u(-)) is a family of compatible 
experiments if there exists (a common) L € £ such that every 7 € l'y, solves (Pr). 


Definition 2.4 (monotonic experiment). An experiment (q(-),u(-)) is monotonic if the control function u(-) is 
smooth monotonic (strictly, 4.e., u/(t) does not vanish for any t). 


From now on we consider monotonic experiments only. In particular, the mapping u + t(u) is well defined, 
continuous and smooth, i.e., belongs to Ü = C©(R;,R). The smoothness assumption, that may look not very 
reasonable, will be justified later (Lem. 3.5). 

In the following, with a little abuse of notations, a time dependent function g is written g(u) instead of g(t(u)) 
whenever necessary. Also, we use the dot * to denote the derivative with respect to time t, and the prime ? to 
denote the derivative with respect to variable u. 


Definition 2.5 (different experiments). Several experiments (q'(-), u?(-)), à = 1,...,4, that visit common con- 
trol values (i.e., such that the ranges of the associated control functions u?(-) overlap) are said to be different 
if there are points # such that g'(t) £ q7(t) for all à Z j, on the overlapping range of control values. 


Remark 2.6. À few obvious but useful observations are in order: 


i) The control of a monotonic experiment can take the value zero at most once. 
P 
ii) Let J be a family of experiments that visit common control values namely, 
HA P ? 


IUwER Vy=(pu()Eel; 314()ER+ | u(t(y)) = wo, 


we may assume without loss of generality that t(7) — 0 for each experiment (shift of time on the 
experiments). 
2.2.8. Statement of the inverse optimal control problem 


The inverse optimal control problem that we shall denote by (x) is the following. 
(aT) Given à family 1’ of experiments, find a cost L in £ such that every 7 € T; solves (Pr). 


3. STUDY OF THE OPTIMAL CONTROL PROBLEM AND ITS INVERSE 


3.1. Study of the optimal control problem (Pr) 


4.1.1. Existence of optimal controls 


We follow exactly the same lines as in [10, Sect. 3.1] for questions dealing with the existence of optimal 
controls. The growth condition L(u) > clu[? with c > 0 and p > 1 used in [10, Sect. 2.1], has been weakened to 
A3 according to [18, Thm. 2.8.1], and thus we can state the following. 

Theorem 3.1. For every go, q1 € Q such that q1 is reachable from qo, there exists a real positive time T* and 
a control function u(-) € L!([0,T*]) such that the corresponding trajectory q(-) solves problem (P;). 
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Notice that Theorem 3.1 does not provide the boundedness of the optimal controls. It will follow from 
Lemma 3.5. 
Remark 3.2. Observe that the existence of a bounded positive T* is not completely trivial and follows from 
the special form of the cost function L. Indeed, it follows easily from assumptions A1 and A3 that there exists 
a positive real number « such that L(u) > @ for every u. Now, if 1 is attainable from qo, let u,(-) : [0,T,] — Q 
be à minimizing sequence. Then, 
aa < C{un(+)) < +00, 


which implies that 7, must be uniformly bounded. 


3.1.2. Application of the Pontryagin Maximum Principle 

Since we do not know yet that optimal controls are bounded in the L® topology, we cannot apply the PMP 
in its classical version. Nonetheless, it is easy to check that problem (Pr) meets all the hypotheses required 
in [18, Thm. 8.7.1] (a more refined version of PMP valid for unbounded controls), and thus we can state the 
following. 
Proposition 3.3. Every solution to (P;) satisfies the PMP. 


In order to apply the PMP to the problem (Pr), we define the following Hamiltonian function: 


R(À,p,q;u) = Ef(q) + Cu + AL(u), 


with p = (6,6) € T5Q and À < 0. 

An extremal curve is by definition a quadruple (X, p(-),g(-),u(-)), meeting the following necessary conditions 
for optimality. 

We recall that the PMP states that if g(-) is an optimal trajectory corresponding to a control u(-) defined on 
an interval [0, T], then there exist an absolutely continuous curve p(-) : [0,7] — T#,Q and À < 0 such that the 
pair (À,p(t)) never vanishes and for à.a. t € [0,T], (q(-),p(-)) satisfies the Hamiltonian system: 


à) = RU PU at) ut) 


(3.1) 
COTON OO) 
and the maximality condition: 
R(A p(6),a(8), u(#)) = max h(A, p(6), a(E), v). (3-2) 


Also, since the final time is free, the Hamiltonian is identically zero along the optimal trajectory, namely, 
E(E)F(a(E)) + AL(u(t)) + C(tju(t) = 0. (33) 


If À Æ 0, the extremal is called normal, while if À = 0, the extremail is called abnormal. Notice that abnor- 
mal extremals are extremal (and often optimal) for any cost since in this case the cost disappears from the 
Hamiltonian. 

Remark 3.4. Here we have to face problems, due to the fact that abnormal trajectories come unavoidably in 
the picture: 


- abnormal extremals, if they are not normal at the same time, cannot provide any information on the system. 
Hence we have to exclude them from our study; 

- abnormals that are normal at the same time are automatically smooth (Lem. 3.5 below). Hence we can assume 
smoothness of our experiments. 

- In fact, and for a purely technical reason, certain trajectories that are not abnormal cause problems: those 
that are abnormal on a piece of themselves only. We shall also exclude them, however, this exclusion can 
certainly be avoided. 
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3.1.3. Description of normal extremals of (Pr) 


In this case, since À 0 we can normalize and assume that À = —1. 
In this case, the maximality condition (3.2) rewrites 


RQ»), a), u(t)) = max k(—1,p(t), alt), v) = EC) F(a(E)) + L'(C(E)), 


with L*(C) = max, (Cu — L(u)) being the Legendre transform of L. The properties of the Legendre transform 
imply the following relations between the optimal control u(-) and C(-): 


u(t) = LY(C(E)),  C(E) = L'(u(#)), (3:4) 
so that the maximality condition (3.3) may be rewritten as 
L(u(t)) — u(t)L'(u(e)) = 6(E) f(a(e)). (3-5) 


The equation for q in (3.1) is nothing but (1.3), while the equation for p, also called the adjoint equation, can 
be explicitly rewritten using € and Ç. Summing up, all the information is contained in the system below: 


= f(x,y) 

ÿ=u= L"(() 

. of 

£— ré (9) (3.6) 


6 = “Lt 
L*(C) + f(x, y) = 0. 


3.1.4. Smoothness of (normal) optimal controls 


Lemma 3.5. The optimal controls u(t) corresponding to normal trajectories of (Pr) are smooth. 


Proof. Assumption A2 implies that L*' is a smooth function. Hence, u is smooth as a component of the solution 
of a smooth differential system. 


Remark 3.6. Here is a place where the strong strict convexity is important. This assumption also implies that 
the set of costs under consideration is an open set. Also, it will be crucial in the normalizations of our costs later, 
and at the end in the possibility to determine uniquely the initial covectors associated with our experiments. 


À by-product of Lemma 3.5 is that in fact our extremal controls are also extremals of the standard (LX) 
PMP (see [1,15]). 
3.1.5. Characterization of abnormal extremals of (Pr) 


As we have already said, certain abnormal extremals will come in the picture. Then we shall characterize 
them. 

First, let us introduce the resolvent R(#, to) of the time-varying linear system (3.6) (third equation) to be the 
matrix solution to À — —R$É, R(to,to) = Ide» (when {6 = 0 we shall write R(t) instead of R(£,0)), and define 
the mapping V through the following lemma. 


Lemma 3.7. To every monotonic smooth trajectory (qo,u(-)), go = (to, Yo), we can associate a well-defined 
smooth map V(-):U —T,,X, defined over the range U oftr ut), such that for every covector £o in R” and 
every u that belongs to the value set of the control, we have 


éoV(u) = E(t(u)) f(a(t(u))) (3-7) 
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Proof. Since the covector é(t) satisfies the linear differential equation £ — —g£ÿf E(t) = éoR(t,to). Ob- 


Dx? 


viously, the matrix RÀ depends only on the experiment so that we can write (£f)(t(u)) = £oV(u), with 
Vu) = R(t(u),to) f(a(t(u)). 
Lemma 3.8. The (monotonic smooth) trajectory t + q{t) = (x(t).y(t)), is the projection of an abnormal 
extremal if and only if span {V(u),u EU} Z TX © R" (or, equivalently, if and only if span {V(u(t)),t € 
(0, 7]} £R”). 

Proof. Assume first that q(-) is the projection of an abnormal extremal. Therefore, and since À — 0, there exists 
a non trivial covector p(-) = (£(-),((-)) such that 


h(0,p(t), (6), u(t)) = 6) f(att)) + C(Eju(t) = max h(0, p(£), a(6),v). 


The maximum being reached in the above equation, we necessarily have C(t) = 0 for all & € [0,T], and then, 
E(t) Æ 0 for every t € [0, T]. Consequently, the zero-hamiltonian condition reads £(t) f(q(t)) = 0, for allt € [0,7] 
or equivalently, for all u € U, 

£(O)V(u) — 0, (38) 


with (0) £ 0. This last equation (3.8) exactly means that {V'(u),u € U} does not span Tr, X. 
Conversely, assume that {V(u),u € U} does not span T,,X. Thus, there exists & € Tr, X such that £oV(u) = 
0 for every u € U, or, equivalently, that for all t € [0, T1], 


é(#)F(a(e)) = 0, (3-9) 


with é(-) being the solution to the Cauchy problem € — —ÉSE(q(t)), £(0) = é0. Differentiating this last equa- 
tion (3.9) with respect to { leads to 


0 = EG) +0 (DEt + EUO ) 


= HN 
Eu, 


which implies that € (-) vanishes identically (by strict monotonicity, the control function u(-) vanishes at most 
once on [0, T|). Hence, the optimal trajectory t + q(t) is the projection onto X x Y of the abnormal extremal t + 
(0,€(6),0, q(t), u(t)). This ends the proof. 


Lemma 3.9. Let y : I — R\ be a smooth curve, defined on some nontrivial interval I. Then, the two following 
assertions are equivalent. 


1. For any subinterval J C I, span{y(t),t € J} = R\, 
2. The set of to such that span {y(to), Y'(to),..., 7 V-1)(t5)} = R\ is (open) dense in I. 


Proof. The direct implication: we prove the result by contraposition. Assume that on some open subinterval 
TCI, span{y(t),..., 7-0 (+); Z R\. For every integer i, define 


r(i) = rank {y(t),... 70 (#)}, 


and let k be the first integer for which r(k) £ k + 1, namely r(k) = k. Let J C Î be a subinterval on which 


span{y(t),...,7%-0(4)} = RF. For every t € J, there exist ao(t), ait), ..., @x-1(t) such that 
k—1 
19 = D autt}0(®. (3-10) 
i=0 
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Notice that the coefficients a;(-)’s are smooth on J. Indeed, relation (3.10) may be rewritten in matrix form as 


ao(t) 
O0 RER P 


ax-1(t) 


which can be smoothly solved in the a;(t)'s (à = 1,...,k — 1) since rank{y(t),.…. 590 00) [te J} = k. Let 
to € J and let R(t,to) be the resolvent of the linear differential equation (3.10), .e., R(t,t0) is the matrix 
solution to the Cauchy problem LR(E, to) = A(t)R(t,to), R(to.to) = Idex, with 


0 1 O0... 0 

A(&)= | 0 DL 
0 ... … 0 1 
ŒD sos se os (071 


We have 
(6,100) = (r(eo), 7/0), 75 0 (0) to)", 


with R(#,t0)* being the transpose of R(t, to). Consequently, for all # € J, 


with Ri(4,to) being the i-th element of the first line of matrix R(#,t0). Therefore, y(J) € R*-! and conse- 
quently y(-) does not span R\ in restriction to J. 

We now prove the converse. If + does not span R\ in restriction to some nontrivial interval J, there exists 
k < N and a k-dimensional subspace EF C R\ such that y(J) € EF. Consequently, for every t € J and every 
integer 4, we have +()(t) € EF. This completes the proof. 


Definition 3.10 (normal experiment). À monotonic experiment (q(-), u(-)), defined on some time interval [0, T 
is normal if q(-) is not the projection of an extremal which is abnormal on [0,71]. 


Definition 3.11 (normal-regular experiment). À monotonic experiment (q(-),u(-)), defined on some time in- 
terval [0, T] is normal-regular if q(-) is not the projection of an extremal which is abnormal on some (nontrivial) 
subinterval of [0, T1]. 


Remark 3.12. Note that a normal experiment maybe non normal-regular. 


By Lemma 3.8, the experiment is normal-regular if V(u) (resp. V(t)) spans R” in restriction to any nontrivial 
subinterval of [uo, ur] (resp. [0, T]). 

As we said, we want to avoid trajectories such that a piece of them is abnormal, (i.e., non normal-regular 
experiments). An immediate consequence of Lemma 3.9 is the following corollary (one has to apply Lem. 3.9 to 
the curve uk Viu)). 


Corollary 3.13. The (monotonic) trajectory, t + q(t) is normal-regular if and only if the set of u such that 
span{Vu),...,V(?-D{u)} = IR is (open) dense in [uo, ur]. 
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3.1.6. Test of abnormality 


A reasonable practical test of abnormality is the standard Gram-matrix test. It may be done either using V(t) 
or Vu). We will write V(r) to emphasize the fact that both can be indifferently considered. Define the Gram 
matrix relative to the interval [71,72] as: 


T2 

Ga =. V(r)V*(r)dr, 
T1 

where V* denotes the transpose of V. 

It is clear that span{V(r),r € J} = RM ïff the symmetric matrix G, is positive definite, iff its smallest 
eigenvalue is strictly positive. Then, the test resumes to decide whether à matrix is positive definite which is 
standard in practice (see e.g. [16] and Ref. [6] therein). 

Differentiating &V(t) as was already done above, we get 0 = £oV'(t) = £oR(t) SE (g(#)). This leads to the 
standard fact that linearization is not controllable along abnormals. 


3.2. Study of the inverse optimal control problem 


3.2.1. Well-posedness of the inverse control problem 


The next lemma will show that the inverse control problem (44) is an ill-posed problem. 


Lemma 3.14. Assume that LE L. Let a, b be arbitrary real numbers with a > 0. Set L(u) = aL(u) + bu. An 
experiment relative to L is also an experiment relative to L. 


Proof. For abnormal experiments the result is obvious. 

Let (q(-),u(-)) be a (normal) experiment relative to L. Then, there exists a covector p(+) = (£(:),6(-)) such 
that the extremal (—1,p(-),q(-),u(-)) is solution to system (3.6). 

Consequently, since L*(C) = sup, (Çu — L) = sup, ((aç + bju — L) = LL*(aç +b), one easily infers that the 
extremal (—1,p(-),q(-),u(-)), with p = (É, ê) = (aë,aç + b), satisfies 


= f(e,u) | 

ÿ=u = L*(() 

£= EST (x,1) 

= FL 
L'(C) + Ef(x, y) = 0, 


showing that (q(-),u(-)) is an experiment relative to L. 


3.2.2. Normalization of the cost function 


Lemma 3.14 shows that the projections onto the base manifold of the extremal trajectories of problems (Pr) 
and (P;r+w) are the same. In particular it says that the cost reconstruction is an ill-posed problem. 

One can only expect to reconstruct a class representative of the cost function under the equivalence relation +, 
L + L if there exists a real a > 0 and a real b such that L(u) = aL(u) + bu. It is worth mentioning that under 
this equivalence relation the minimum of the cost function is not preserved but the set of control values such 
that L(u) —uL'(u) = 0 is preserved (indeed, due to assumption A1, A2, A3, there exist exactly two values for 
which this relation holds: draw the tangents to the graph of L through the origin). 

The following lemma gives a way to select a special class representative for the cost function. 


Lemma 3.15 (resetting lemma). At uo the cost function L to be reconstructed can be normalized 80 
that: L'(uo) = 0 and L''(uo) = 1. 
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Proof. Let L be a class representative of the cost we aim to reconstruct. In order to get the required normal- 
ization, we need to solve for à, b real numbers, with a > 0 the following system 


L'{uo) 1\ fa\ _ {0 

L' (uo) 0 b}  \1/° 
which is always possible for a cost function L in L: the coefficient a obtained in this way is positive (a = 1/L"(uo)) 
due to the strong convexity hypothesis. 


Remark 3.16. Notice that our normalization does not preserve the set £: after normalization, the assump- 
tion A1 may not hold, but this will be of no importance in the following. In fact the right way to define the 
class £ would be to replace the positiveness of L by the weaker L(0) > 0. But this is unnatural for the purpose 
of proving the existence of minimizers (Thm. 3.1). 


3.2.3. Eristence and uniqueness of the reconstructed cost 


Lemma 3.17. Fix a system f in F and a smooth admissible trajectory (q(-),u(-)), with u(-) monotonic and 
defined on some time interval [t1,t2]. Fix a covector £o € R"\{0}. Then, there exists a unique smooth function L 
defined over the range oft - u(t), and such that L'(u(to)) = 0 for some arbitrary point to € [t1,t2], for which the 
pair (q(-),u(-)) is the projection of a normal extremal of problem (Pr) associated with the covector (£(-),6(:-)), 


(£(ta), C(ta)) = (60,0). 


Proof. Choose to € [t1,t2] and set uo = u(to), go = g(to). Let V(-) be the function associated to (g(-),u(-)) and 
defined over the range of u(-) accordingly to Lemma 3.7. 


We begin with the existence of L and we also give a formula for it. The range of u is [u1, u2]. Assume wo # 0. 
Consider any a € [u1,u2]. The reader can easily check that 


L(u) = ÉoV(uo) + (éoV (wo) — éoV(a)) (2 - 1) — / j COOP (3.11) 
is such that L'(uo) = 0, and satisfes 
Lu) — uL'(u) = &V(u). (3.12) 


If 0 € fui, u2], this L(u) is smooth and is the unique solution of the standard ODE (3.12) with L(uo) = 
£oV(uo), L'(uo) = 0. If 0 € [u1, u2], let us chose a = 0. Define the function g by g(u) = £oV(0) — éoV(u). We 
have 


g'(u) = OO AO 


Hence, g(0) = g/(0) = 0, and we may write, g(u) = u?g(u), with ÿ smooth. L(u) — £oV(uo) + (£oV(uo) — 
PTS = LT us, g(v)dv is a smooth solution of (3.12), with again L(uo) = &oV(uo), L'(uo) = 0. 
Suppose now that uo = 0. Then, set L(u) = £oV(uo) + Die foV(0)—éoV(e) dv. Again, L(u) —uL'(u) = £oV(u), 
g(u) = w?g(u), with ÿ smooth. L(u) = £oV (uo) + g(v)dv, is a smooth function that meets L(0) = £6V(uo), 
L'(0) = 0. 
Note that the solution L(u) is well defined and smooth over the full interval [u1, ul]. 
We now prove uniqueness, when 0 € [u1,u2]. We assume that 0 € [uo, u2] (the case 0 € [u1, uo] is similar). 
Let Li and Z2 be two solutions to equation (3.12) such that L'(uo) = Li(uo) = 0. Then, # = Li — La is 
a C'® solution to 


vu) En uy'(u = 0, y(uo) = 0, Y'(uo) = 0. 


If wo is nonzero, then 4 (a smooth function) is identically zero on [uo, 0]. On ]0, u2], we must have Y(u) = Ku, 
and since we look for a smooth function, d'(u) — K must be zero. At the end, 4 — 0, which proves uniqueness. 
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It remains to prove that the pair (q(-), u(-)) is the projection of an extremal trajectory of the problem (P;). 
By hypothesis, the given trajectory is an admissible trajectory of system (1.3). 

By definition, V(u(t)) = R(t,to)f(a(t)). Set £(t) = ÉoR(t,to) so that é(-) satisfies the third equation of 
system (3.6). 

Set C(t) = L'(u(t)). Thus, C(to) = 0, and 


6(E) = L'(u(t))à(®) 


eu) DE (a(t()) (D 


- (4) SE (alta) 


which shows that ((-) satisfies the fourth equation of system (3.6). 
By construction equation (3.5) is satisfied. This ends the proof. 


Remark 3.18. Note that it could be that the trajectory under consideration in the lemma is abnormal. In that 
case, by (3.11), the smooth ZL reconstructed in the lemma is a constant. 


3.3. Validation of the hypotheses 


3.8.1. Validation of the hypothesis L(u) 


We strongly think that the best way to numerically validate the hypothesis L — L{u) is to solve the set of 
equations (4.3) in the least squares sense. The interest of the least squares approach is that it allows an arbitrary 
number of experiments (strictly larger than 2). The least squares solution provides adjoint initial vectors, that 
can be used to reconstruct the cost over the set of values of u visited during all experiments (provided that they 
overlap). 


Then as soon as the experiments are compatible (4.e., the trajectories are actually solutions of an optimal 
control problem of the required form), the least squares procedure provides a solution. Moreover, if the answer 
is positive, we have a constructive way to reconstruct the cost. 


3.8.2. Testing the strong convexity of the cost 


Once the covector £ has been found, it is easy to check the strong convexity of the cost function. Indeed, 
differentiating equation (3.5) with respect to w leads to 


uL'"(u) + £V'(u) = 0. 


Thus, when u £ 0, the cost function is strongly convex at u if 6 V'(u) < 0. At u = 0, it is enough to differentiate 
equation (3.5) once more to see that the cost function is strongly convex at zero if 60V” (0) < 0. 
Consequently, if the system of equations (4.3) admits a solution, it is reasonable to add one of the two above 


inequalities in order to validate the strong convexity of the cost. 
4. THE RESULTS 


4.1. Statement of the main theorem 


The inverse control problem being properly posed, we can state our main theoretical result. To state it, let 
us endow the set F of systems under consideration with the C’©7 Whitney topology. 
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Theorem 4.1. There is a residual set of systems in F such that three monotonic, different, compatible and 
normal-regular experiments which visit common control values are enough to reconstruct the cost function L 
over the set of values of uw visited within the three experiments. In other terms, cost reconstruction with 
three experiments (whose range of controls overlap) is a generic property. Moreover, generically, as shown 
in Theorem 4.15, two experiments may not be enough. 


Remark 4.2. Of course, it does not mean that, in practice experiments are monotonic. À non monotonic 
experiment consists of several monotonic pieces. For instance, an experiment which visits three times the same 
value of the control is generally sufficient for the cost reconstruction. 


The proof of Theorem 4.1 is postponed to Section 4.3, after certain preliminaries given in the next section. 


4.2. Preliminaries for the proof of Theorem 4.1 


4.2.1. Notations 


We will use the following notations in the remaining of the paper: 

Qc = {ao =(q',...,a)1at Ag for 1<a<B<{}. 

e R*!T1 denotes the set of typical {-tuples (4},...,Uf) of k-jets of controls, Hi — (ui, u} = ü’(0), ..., 
uÿ, = uw") (0)), meeting uÿ = uo £ 0, u}Z 0,... ,uf# 0 (common nonzero value wQ of the control at t — 0, 
and first derivatives of the control nonvanishing at t = 0). 

e Zke = Q(o X R£‘+1, Hence, dim Z4,e = ln + 1 + k£. Independently of the values of £ and £, an element of 
Zr.s is denoted by z. 

e The bundle (EF)! is the restriction of the product bundle (J*F)* to Qt, where JF denotes the bundle 
over Q of k-jets of elements of F. Typical elements of (JÉF) are tuples 


Ffotawo) = Ghf..., jf), qw =(a",...,4) € Q(. 


4.2.2. Estimation of the number of needed experiments for the cost reconstruction 

For fe F, let (q!(-),u!(-)),...,(q{(-),ut(-)) be £ (£ > 1) compatible, monotonie, different and normal-regular 
experiments to which correspond initial covectors €, .…. ES respectively. Denote by U!,...,U‘ the range values 
of ul(-),...,ut(-) respectively. 

Assume that these experiments visit a common control value wo; and thus, by monotonicity, visit a common 
control open interval U C R. Without loss of generality, we may assume that 


uo À 0; 

q*(0) £ qg°(0) for 1 < a < B £ L (since the experiments are different); 

üuo = ul(0) = --: = ut(0) (up to a translation of time on each experiment); 

üù(0) Z0,...,ù(0) Z 0. 

All the matrices A? = (V® (ao), V%(uo),..…, va (uo)) have full rank, where V®(.) is the map defined in 


Lemma 3.7. 


Remark 4.3. - The reason why we assume that wo # 0 will be made clear later, in the proof of Lemma 4.12. 
It is related with the fact that, when wo = 0 (which is possible), in formula (4.10), Lemma 4.12, £V®(uo) 
appears multiplied by wo. 

- The fourth assumption is justified by (strong) monotonicity of the experiments. 

- The fifth assumption is a consequence of the fact that the experiments are normal-regular and of 
Corollary 3.13. 

- The reason for the second among the five preceding items, will be made clear later, in the proof of 
Corollary 4.13. It is related to the fact that we will need to fix k-jets of the same function at different points 
simultaneously. 
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On the interval U of common control values, the following holds: 
&V'(u) =: = EV'(u) = L(u) — uL'(u), (41) 
and the non homogeneous equation 
66" (uo) = --: = 66V* (uo) = —u0 (0), (4.2) 


has to be considered if we wish to reconstruct the cost function normalized according to Lemma 3.15. 
Equations (4.1) and (4.2) imply that for all integers k > 0 


ss (= 0, 43 
V1 (u0) = —uo, 
where M(2) is the {n x ({ — 1)(k+ 1) matrix defined by 
A7, (Ù à 
—A? A? 
M(2)= | 9 A3". Oo |: (4.4) 
. 4 A 
0 ... O0 —A4 


with A? = (VS (u0), V#/(uo), …, Va) (u0)) being the k-jet of V®(-) at wo, and z has been defined in Sec- 


tion 4.2.1 above. 


Remark 4.4. Notice that the equations &V%"(uo) = —uo for à = 2,...,{ are omitted since they follow from 
the second equation of system (4.3). 


The main theorem (Thm. 4.1) is a consequence of the fact that system (4.3) has (gencrically) a unique 
solution, as we shall prove. 

Before, let us examine the different cases showing up for the solutions of system (4.3). Let Sy ; denote the space 
of Un x (£ — 1)(k + 1) matrices of the form (4.4). Assume that {n < (4 — 1)(k +1). This is not a restriction since, 
later, k will be taken large. The set of matrices Sxs is stratified by the rank, each strata being a submanifold 
of the set of {n x ({ —1)(k+ 1) matrices. Namely, 


ln 


She = U Sretr), 


r=0 
with Sy e(r) being the subset of matrices of corank r. 


The three following situations may appear: 


M(z) € Sx,e(0). In this case M(z) has full rank. The unique solution to the first equation of system (4.3) is 
zero and consequently the system (4.4) has no solution. This may happen if: either the considered experiments 
are not compatible (they are not extremal for the same cost), or the assumptions made on the cost function 
are not satisfied. 


Remark 4.5. This case cannot happen under the assumption that the experiments are monotonic, different 
and compatible. 
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M(z2) € Sxe(1). This case is the one we are interested with. We will show immediately that system (4.3) 
admits à unique solution and then, according to Lemma 3.17 the cost is uniquely determined. 
Let us start with &o = (£,...,€f) Z O vanishing on M(2). If EV°/(uo) Z 0, the inhomogeneous condition 
in (4.3) can be satisfied just multiplying & by an appropriate nonzero constant, and the solution so obtained 
to equations (4.3) is unique. 
If ES V%(uo) = 0, let us show that we can change wo (an arbitrarily small change for the 5 open conditions 
of the beginning of the section still hold true) for ES V®%’(uo) Z 0. 
The experiment being assumed normal-regular, we claim that € V®%’(u(t)) Z 0 for some t in an open interval 
arbitrarily close to t(uo) = 0. Indeed, assume that &°(t) = ES R°*(t) f(q*(t)) is constant on some interval 
(the resolvent R*(t) has been defined at the beginning of Sect. 3.1.5). 
Then, if this constant is zero, it means exactly that the trajectory q*(-) is abnormal on this neighborhood, 
which contradicts the normal-regular assumption. If this constant (say k) is non zero, it means that the 
trajectory is normal and 

EV®(u) = L(u) —uL/{u) =k #0, 

with, in addition, L'(uo) = 0 and wo # 0. It yields (integrating the differential equation) to L(u) = k 
contradicting the strict convexity of L. 
Then, we can take for uo a u*(fo) such that &*(40) £ 0. The new £5 will be automatically (a scalar multiple 
of) é&R°(t0). Moreover, M{z(t)) will still be in Sx,s(1): it cannot belong to S; (0) since the experiments 
are compatible, and $4s is stratified by the rank. 

M(2) E U,z2 Sre(r). This case is the one we want to avoid. Indeed, if corank M(2) > 2, then system (4.3) 
admits a solution but this solution is not unique and, according to Theorem 4.15, the cost L is not uniquely 
determined. 


In the following we will study this last case and prove (using Thom’s transversality theory) that it is generically 
impossible as soon as { > 3 and k is chosen large enough. 


4.2.8. The bad sets 
Let M4 be the set of structured matrices of the form (4.4). 


Definition 4.6. Let Bxv(r) be the subset of My. which elements M satisfy: 
corank M = r > 2; 
every submatrix AŸ has full rank (a = 1,...,4). 
Set Bx.e — LÉ Brilr). 
Lemma 4.7. If M € Bye, its cokernel has no non zero element of the form (0,£2,...,€t) or (£1,0,€%,...,€t) 
un, 20m (E,22.060 5 0, 


Proof. We prove the result by contradiction. Assume without loss of generality that (0,£2,...,€*) € coker M. 
This means that (0,€2,...,£t)M = 0. Consequently, 2 4? = 0, which implies that £? = 0 since 4? as full rank. 
By induction, we get £* = 0 for all «a = 1,...,4. 


Let Gr(r,R°) denote the Grassmannian of R°”, that is the space of all r-dimensional vector subspaces of R!”. 
With a little abuse a notation, we identify a r-dimensional linear space I € Gr(r,R°*) with any r-tuple vectors 


(É1,..., 7) = (EE, 5er) euch that IT = span(éi,..….,ér). 


Definition 4.8. Let Bru(r) be the subset of Gr(r, R°*) x M4.e of all tuples (&1, at É, M) such that M € Bx4 
and (£1,...,&)M = 0. 


Notice that the set Bx4 is the projection of Be = CES Byu(r) parallel to the Grassmannian Gr(r, R{?). 
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4.2.4. Estimation of the codimension of By + 


The Thom’s coranks formula (see [17]) does not hold in general for the structured matrices M (1.e., of the 
form (4.4)), but it still holds for the matrices M] whose blocks have full rank. More precisely we have the 
following lemma. 


Lemma 4.9. The semi-algebraic set Bx4 satisfies 


codim Bys(r) >r((L—1)(k+1)—-{{(n—r)), (4.5) 
where we have assumed that U(n — r) < (£—1)(k+1). 


Proof. Let M € Bye(r). Then, we can find r independent vectors Ë = reset Ë, = (£l,...,ét) in the 
cokernel of M, 1.e., 
ES AN ETAT LG, ÿ=1,...,r a—=1,...,£4—1, 


where, for simplicity, we have set A = AŸ (a = 1,...,4). We may assume without loss of generality that, 
for a = 1,...,{ —1, the &%°s (à = 1,...,r) are the rth first vectors of the canonical dual basis of R”, namely, 
that €® = (0,...,0,1,0,...,0). Indeed, if €,...,£l were dependent, there would exist (A1,...,1;) Z (0,...,0) 
so that Al +... + À,€l = 0. Consequently, we would have a non zero element of the form (0,£2,...,£,) in the 


cokernel of M contradicting that M belongs to Bx(r) (and the same holds for Éa, a=92,...,1-— 1): 
Now, in a neighborhood of (£,...,£&, M), let us consider the set of r({ — 1)(k + 1) equations 


F9 = (E8 + 668) (AT + SAT) — (EXT + ER TT) (AL 4 GATE) LO, (4.6) 


2 


(with à = 1,...,r and à = 1,...,{ — 1), which we want to solve in w = (u1,w2,...,W,(g_1)(&+1)) = 
(8Af,...,0A%) where 6 A denotes the i-th line of 5A°. 
Notice that, 


0F* Ô Ô Ô 
à ne. Fe nn a A — a+1 A+ = 
Le 1 CE) _ SFA —5et15 


AË = GAS 
Ow ô 2 Ô 2 


=0 w w=0 


implies that the Jacobian matrix at zero with respect to the variables w of the mapping F = Fi x Fi x... x 
F1 x F? x... x Ft associated to the system (4.6) has the form 


Idgx+1 * se * 
OF 0 Idgx+1 É 
— (0) = 
Du ) 
D à Cd 


Hence, by the implicit function theorem, we can smoothly solve system (4.6) in a neighborhood of zero. 
Therefore the set of solutions to system (4.6) is a manifold of codimension r(£ — 1)(k + 1), namely, 


codim By.o(r) = r(£ — 1)(k +1). 
Since Bz.(r) is the projection of Bx4(r) on My, parallel to the Grassmannian Gr(r,R{"), we have: 


codim By.s(r) — codim Gr(r,R°*) 
r(£—1)(k+1)—r(£n—r). 


codim Bye(r) 


This is the Thom’s “product of coranks” bound. 
Since r > 2, we get (k large): 


codim Bz,e > 2{((L — 1)(k +1) — Un +2). 
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The following proposition, which is sufficient to prove our main theorem (Thm. 4.1), is à consequence of the 
considerations of this section. 


Proposition 4.10. If { = 3 and k is large enough, then the set of f € F such that the map z -— M{(z) avoids 
B%,3, in restriction to 7% 3, is residual in the Whitney topology. 


Remark 4.11. The estimate codim B% 4 © 2(4 — 1)k for k large shows that we cannot expect that two experi- 
ments are enough: see Sections 4.3.3 and 4.3.5. 


4.8. Technical mathematical tools and final proof 


In this last section, besides à technical lemma and the final proof of our main theorem, we show that in 
general two experiments are not enough. 


4.8.1. À crucial lemma and its corollary 


Here, R*+1 denotes the set of (k+1)-tuples (uo,u1,...,ux) with wo, u1 both nonzero (k-jets of smooth controls 
that are nonzero and monotonic). Define the following mapping: 


SAP RTL RIRE TN) 
°F j0u) + ao V 
Lemma 4.12. The mapping = is à surjective submersion. 


Proof. Note. Along this proof, we will make a constant abuse, for simplicity in the notations: we are treating 
k-jets of systems f and controls u(:), as the system f or the control u(-) themselves. This is justified by the 
fact that all the manipulations we do depend only on the k-jets of the objects. For instance, the k-jets of the 
function t(u) depend (smoothly and diffeomorphically) on the k-jets of u(t) only, due to the strict monotonicity, 
reflected in the definition of R£+1. 


A point (uo,u1...,ux) € R**l is given, together with a k-jet j"f € J*F with source go = (0, yo). Then, 
(V(uo), V'(uo), .…,V ®(uo)) is well defined. 

Assume that q = q(t) (equivalently q(u)) is a (monotonic smooth) trajectory defined on a neighborhood of 
zero (equivalently wo) such that q(0) = qo (equivalently q{(uo) = go). Let R(u) = R(t(u)) be the resolvent defined 
at the beginning of Section 3.1.5. 

Below we shall compute the successive derivatives with respect to u of V(u). Namely, 


V'(u) = R(u)rx. (ut/(u)adr Fo(g)). 
Vu) = R(u)rx ((ut/(u))" ad, Fo(g) — u (#/(u))" adr, Fi (9) 


—(r1) (0) l'adfetFi(g) + ub.(u) + di(u)), 22, (47 


with mx. being the tangent map to the canonical projection mx : X x Ÿ — X parallel to Y, and where the 
notation ,(u) means à term linear in Lie brackets (evaluated at q) of length not greater than r+1 containing F6 
and F; and in which F; and F; appear at least once (r > 1) and at most r — 1 times (r > 2). The coefficients 
of the linear expression #,(u) are functions that belong to R[u, t'(u),...,#("){u)]. 

Let us now prove formula (4.7). Introduce the following notation for Lie brackets 


Fr =[F,,[F,,...[8,.,,F,]...]] T = (i1,...,%), Hi=iü +... +i. 
Define the set Z,. by 
r+1 
D = U {re 10,1} | 1< | <r-1, 1<t-l<r-1}. 
£=1 
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Then, 
Vr(u) = D crr(u)F1(g) ETX X {Ory}, cru) € Ru,t'(u),.…., #0) (u)]. 
IET, 
For any 0 € T5, X, let pt) = (é(),6(-)) be the solution to 5 — —p£r — upiE , P(0) = (éo, 0). 


Consequently, taking into account that cr, (u)' is of the form er,;+1(u) (which is obvious), and that w,(u) € 
TX x {0ry}, we get 


érxl(u) = p#.(u) (48) 
= SE és) F(o) + crr(u(n) © (p Fi(a) 
TEL 
_ > Cr (u)p Fr (q) Li cr,r(u)t'(u)p adpo+uF Fr (g) 
IET; 
TC) (4.9) 


We now prove formula (4.7) recursively. We have 
éoV'(u) = (6rx.Fo(a)) (u) 
(p Fo(g)) (u) 
d / 
= + RG) (u) 


= padr,+ur Fo(g)t'(u) 
= érx.adr, Fo(g)ut'(u), (4.10) 


which gives the formula for r = 1. 
Differentiation of formula (4.10) with respect to w gives the formula for r = 2: 


&oV"{u) = (padr Fo(q)ut'(u))' (u) 

= À (padr, Fa) d'ut/(u) + pad Fo(a) (ut/(u)) 
Erx+ (aa rt ad, Fo(q)t(u)ut/(u) + adr, Fo(q) (ut'(u))") 
= Enx ((ut/(u)Ÿ adé, Fo(g) — u (#(u)) ad? F1(9) — d'(u)adr, F1(a) + uvalu) + 41(u)), 


with dou) = t’(u)adr, Fo(q), and Yi(u) = 0. 
Now, differentiation of formula (4.7) with respect to u and consideration of equation (4.9), leads straightfor- 
wardly to formula (4.7) at rank r +1. 
Note that since F5 = (f,0), F1 = (0,1), we can rewrite the formula (4.7) at u = wo as: 
d'V 


(ue) = (uot/(ue))" D (a) + 6 (uo) 


where @-(uo) depends on successive derivatives of f with respect to y up to order strictly lower than r. And 
since 4o, t’(uo) are both nonzero, the result follows. 


! 


Corollary 4.13. The mapping 
Et: (EF x RÉ = Mie 
(fo... Uf) + M(2), 
with z = (q!,...,q!,U},...,UE), is a surjective submersion. 


Proof. The proof is an immediate consequence of Lemma 4.12, since for every q(4) € Q(), the L points Ge Q 
are distinct. 
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4.8.2. Proof of Theorem 4.1 


In this section, we will apply a slight modification of the standard multijet transversality theorem (see e.g. [12, 
Thm. 4.13]). This theorem is stated for bundles of jets of smooth mappings, but it holds also for jet bundles of 
sections of a bundle. This is known, and completely clear, since the arguments are essentially local, and locally, 
a section of a bundle is just a smooth mapping to the fiber. 

We state the theorem in our own context only: 


Theorem 4.14. Let W be a stratified subset of (JFF)£. Let Tw = {f € F | j* fe) À W}. Then, Tw is a 
residual subset of F for the Whitney topology. 


4.3.8. Definition of Nye 
The set Ne C JF that we will construct below, will, in some sense, represent the set B% in IF 


By Corollary 4.13, the map ©° is transversal to the points. Therefore, (£*)-!(B,,) is a semi-algebraic 

Whitney-b stratified set of the same codimension cg, t.e., cx,.e > 2((0 — 1)(k + 1) — {n +2) by Lemma 4.9. 
Let Ti : (IF x RétTI — (FF)! denotes the canonical projection on the first factor parallel to RÂ!+. 
Define the set M4 4 by: 


Nre = m (5°) (Buse) 


Consequently, N,.4 is a stratified set of codimension: 


codim By,e — (kl +1) 
2((4—1)(k+1) —-£n +2) —Kk£—1 


codimNye > 
2 


For £ = 3, 
codimNy3 > k —6n +7. 
4.8.4. End of the proof of Theorem 4.1 


By Theorem 4.14, the set G of f € F such that TE) is transversal to N43 is residual for the Whitney 
topology. If k is large enough, it means that jf) avoids Ny 3. 

But avoiding Nx3 means exactly that the map z + M(z) avoids B43, in restriction to Zx3. By Proposi- 
tion 4.10, this ends the proof of Theorem 4.1 


4.8.5. Necessity of three experiments 
Theorem 4.15. Generically, two experiments are not enough in order to reconstruct the cost. 
Proof. Let (q4(-)},u4(-)) and (gÆ(-),u?(-)) be two compatible, monotonie, different and nee ex- 


periments to which correspond V4 and VF. According to Lemma 3.9, there exists uo Z 0 such that j, se k y/A 
and 1 VÆ have full rank. By compatibility of the experiments, there exist £À and £? such that 


A A 1 (4.11) 
EVA (uo) = —uo. (4.12) 
Set M(z) = (j5,V4, —j5 VÆ)* and suppose that dim coker M(z) = 2 (the strata of matrices M(z) whose 


cokernels have dimension strictly larger than two is generically avoided, by formula (4.5), but the stratum 
corresponding to corank 2 is not avoided in general). Then, the set Æ = {(£4,£F) | (4.11) and (4.12) are 
satisfied} has dimension one. Hence, there exist two independant pairs of covectors (££},£P) and (£4,€4) such 


that E — Gr.) + R(E. 67). 
By (4.12), eA > are all nonzero. Assume that £A,€f\ are linearly dependant, then, some Éa D 


Set Ép = £P + X£P, and Ê = (0, 6n) € E for êg £ 0. This again is impossible by normal-regularity, since it 
implies éB JV = = 0. Therefore, €, é! are linearly independant. 
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The cost function defined by 


u qrA _vA( 
LA) = (E4 + XEA) (v4{uo) tu f TT O4) 


0 
is clearly a solution of the inverse problem for any À in R. PE L(À,u) depends on À unless ei 

But, L(\,u) — uPEQuu) 1) = = (£f + XES)V A(u), hence if $E = 0 then differentiating w.r.t. À, 0 = Z(L(\u) — 
PEER DL) £$V A(u), which again cannot be identically zero by normal regularity. 
Condeaentie two experiments do not allow to reconstruct the cost function uniquely. 


u 


5. THE INVERSE OPTIMAL CONTROL PROBLEM FOR THE DUBINS SYSTEM 


5.1. Preliminaries 


This part is concerned with the resolution of problem (44) for Unmanned Aerial Vehicles (UAVS). 
In this study, we consider a HALE (High Altitude Long Endurance) UAV flying at constant speed and 
altitude. The equations of motion are the ones of the Dubins car (see [9] for a justification), .e., 


à = cos Ü 
ÿ = sin 0 (5.1) 
ô = u, 


with q = (x,y,0) € R? x S! being the state (where (x,y) € R? is the UAV's coordinate in the constant altitude 
plane, and 0 the yaw angle), and u € R being the control variable. 

Inspired from the works [6,8,10] we assume that the trajectories flown by experimented pilots are solutions of 
some optimal control problem. The problem is to decide what is minimized. As written in [8], this assumption 
is based on a “nowadays widely accepted paradigm in neurophysiology which says that, among all possible 
movements, the accomplished ones satisfy suitable optimality criteria (see [14] for a review)”. In [8,10] the 
authors address the problem of reconstruction of the cost minimized in human locomotion. It's worth to notice 
that our problem is very similar since HALE drones behave kinematically more or less as a human being moving 
on à plane (constant altitude, constant speed). 


Remark 5.1. 1. Up to a suitable rototranslation and a time shift, we may assume that optimal trajectories 
meet g(0) = (0,0,0) (invariance under motions). 

2. Problem (Pr) for the Dubin’s model does not admit any abnormal extremal as the reader can easily check. 

3. Note that u is the curvature of the trajectory in the (x,y) plane. A special class of trajectories consists 
of those trajectories with a single inflexion point (1.e., a zero curvature point). Although our method allows to 
care about all trajectories, these are of particular interest, as will be discussed in Section 5.3. In that case, it 
would be completely natural to consider that wo — 0, as in the papers [8,10]. In fact, uo — 0 is a case that we 
tried to avoid at several places in this paper (we overcome this case by small perturbation). Therefore, we will 
not make this normalization here. 


5.2. Main results in Dubins case 


We seek for a cost function normalized as in Lemma 3.15. For system (3.1), the free time condition writes: 


L(u(t)) — upo(t) = pr(t) cos O(t) + py(t) sin 4(t), (5.2) 
and the adjoint equations are: 
Dy = 0 (5.3) 


Do = Pr Sin 0 — py cos 0. 


Therefore in the remaining of this section, p4 and p, are real constants. 


123 


20 ALAIN AJAMI ET AL. 


Taking into account equation (3.4) and the normalization of Lemma 3.15, the evaluation at { = 0 of the last 
equation of system (5.3) gives 


(0) = L’(u(0))à(0) = pa(0) = —py. (5.4) 


And the last equation of system (5.3) may be rewritten 


DO = Paÿ — Pyb. (5.5) 


Since æ(0) = y(0) = pa(0) = 0, we get 
po(t) = pay(t) — pyx(t), 
which, according to (3.4) and (5.2), leads to 


L(u(t)) = px (cos 0(t) + u(t)y(t)) — (0) (sin 0(t) — u(t)x(t)). (5.6) 
In particular, the cost L is known as soon as the constant p, is known. 


5.2.1. Consider a single experiment O(u) 


Let us consider a single monotonic experiment #(u) = #(t(u)) parametrized by « (the curvature). The choice 
of any arbitrary p, determines L(u) according to equation (5.6). Since the L(u) so obtained depends on p,, we 
conclude that one experiment is not enough. 


5.2.2. Consider two compatible experiments O(u), O(u) 


Assume that the respective domains U and Ü of these experiments have a non void intersection. It is always 
the case for (x, y) trajectories with an inflexion, to which we can apply the previous normalization g(0) = (0, 0,0), 
in particular 4(uo) = O(uo) = 0, which is responsible for (5.4). 

Since the two experiments are compatible we have for all u in U ( Ü 


L(u) — uL/(u) = p; cos O(u) + p, sin O(u) 
= px cos O(u) + py sin O(u). 


Therefore, we must have for all u in U 
px cos O(u) — px cos O(u) = —p, sin O(u) + p, sin O(u). (5.7) 
We have the following lemma. 


Lemma 5.2. For two (strictly) monotonic compatible experiments, the set of equations (5.7) determines px 
and p; as soon as the functions cos(6(-)) and + cos(0(-)) do not coincide on their common domain. 


Proof. For two arbitrary control values u1, u2 in U, define the determinant: 


D(u, u2) = det (Cu ) | ( . ee ) | 


Let us assume that D(u:,u2) = 0 for all u1, u2 in U (\ Ü. We deduce that, since 0 and Ô play the same role, 


2a€R, Vu eu cos Ü(u) = a cos O(u). 


By substitution of this last equation in (5.7) we get, for all u in U ( Le 


Px cos O(u) — pra cos Ou) = —pyey 1 — cos? Ou) + pyév 1 — a? cos? Ou). 
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with e,ë € {—1,1}. This last equation implies that 

P(X)=AX?-92BX1C—0, (5.8) 
with 
— cos? O(u), 
…. - 2 2- \2 = 40 2- \2 

(Pr — Dx@) + (Py — @ y) (Px . Pxa) LS (Py +4 Dy) ; 
3 - : nt 

(5 — 55) (ps — d°5,) + (p} + 55) (pa — Para), 

2 _ -2)2 
(p, —5,) : 

Note that p, À 0 by (5.4) and strict monotonicity. If one among the coefficients À, B, C' is different from zero, 
this implies that X is a constant, £.e., cos O(u) is a constant over UA\U. Therefore O(u) and @(t) are constant 
over some nontrivial interval, this is impossible by strict monotonicity. 

If À, B,C are all zero, this implies p, = +ÿ, by C = 0. Plugging in B = 0 gives p+ — pra = 0, since py # 0. 


Plugging both in À = 0 gives a = +1. : 
Hence, either for some w1, u2 in UAU, D(w,u2) # 0, and p;,p4 are uniquely determined or cos(#(u)) and 


+ cos(b(u)) coincide on U AU. 


Q & à» x 
| 


By the previous lemma, given two monotonic compatible experiments (renormalized as above, 1.e., such that 
q(0) = (0,0,0)), there are two cases: 


1. the constants p,,p, are uniquely determined, or; L 
2. cos(O(u)) —+cos(O(u)) on some nontrivial interval. This implies that O(u) —+O(u) (the situation O(u) — 
x + O(u) does not happen since 4(0) = 4(0) = 0). 


Case 1 implies that L(-) is uniquely determined by (5.6). 

Case 2 implies that, after normalization g(0) = (0,0,0) of the experiments, u(t) = &ä(+t). This is equivalent 
to the fact that both experiments are either the same or deduced one from the other by an isometry of the 
plane, that changes orientation (see Rem. 5.5 below). 


Definition 5.3. We say that two experiments are M-different, if their x, y trajectories cannot be deduced one 
from the other by an isometry of the plane. 


We have proven the following theorem: 


Theorem 5.4. Considering system (5.1), the problem (44) can be solved if the set of experiments is composed 
by two monotonic, M-different and compatible experiments which visit common values of the control. 


Remark 5.5. One can read, in (5.6) that if u(t) is changed for u(—t),x(t) is changed for —x(—t), O(t) is 
changed for —40(—t), then L(u) remains unchanged. It means that we cannot expect to determine L(-) from two 
experiments deduced one from the other by a reflexion w.r.t. some axis. This corresponds also to the trajectory 
(&(t),g(t)) being the trajectory (x(t), y(t)) followed in reversed time. 


5.3. Numerical results 


The theory exposed in the preceding subsections was tested by means of both simulated and experimental 
data. More detailed simulation results will be presented in a forthcoming paper, relative to our global project 
on drones. Roughly speaking, the reconstruction algorithm is, as we explained in Section 3.3, based upon least- 
squares applied to equation (5.6), in order to reconstruct the adjoint vectors. Once the adjoint vectors are 
known, the cost can be reconstructed on the full domain visited by all experiments. 
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Cost reconstruction 


#  Theorical cost 


FIGURE 1. Ideal reconstruction 
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FIGURE 2. Reconstruction (left) on the basis of two experiments from a pilot (right). 


Figure 1 shows just reconstruction based upon perfect experiments, computed from the known cost. Of course, 
this is just checking that the theory is not false, and that no purely numerical problem appears. This last point 
is not so obvious, as explained just below. 

On the other hand, experimental data were obtained from the simulator (developed in our project) [2] in its 
manual mode. Datas shown here, in Figure 2 come from a beginner pilot. The pilot was just asked twice to 
change from flyving direction for à parallel one, which implies the existence of an inflexion point (zero curvature). 
For the sake of testing the algorithm performances, we considered minimum possible level of information, 1.e., 
two trajectories only. Both trajectories were restricted to the maximum monotonicity interval. 

One can see that, in this type of experiment (change for a parallel way), monotonicity looks not satisfied (it is 
very close not to be, and here, there is a real sensitivity problem in practice). Two costs have been reconstructed 
from the determination of both adjoint vectors, obtained from least squares. 

One can see that this pilot actually looks to minimize a cost of the requested form, and this cost looks 
properly reconstructed (both reconstructions are very close). 

Figure 3 shows a reconstruction, from the same pilot, on the basis on two trajectories without inflexion points. 
These second experiments have been performed one month later. 
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FIGURE 4. Comparison of both reconstructed criteria. 


In Figure 4, we show a comparison of the reconstructed criteria for the two previous couples of experiments. 
Of course we make the best possible fitting using our correction term L(u) = a(L(u) + bu). It seems that this 
pilot really minimized the same cost. 

As a conclusion, the first trials we made seem to confirm the applicability of the method, and the reliability 
of the algorithm. 

À measurement series with experimented pilots is on the point to be performed on their own training simu- 
lator. 


6. EXTENSION OF THE RESULTS TO THE CASE q = F(q) + uF(q) 


In this section, we generalize the results obtained in the previous sections to the more general class of smooth 


systems of the form 
9 = Fo(q) +uF(q), qe Q, uER, (6.1) 


with the state space Q being an n-dimensional connected, Hausdorff and paracompact smooth manifold. 
To do so we follow exactly the same lines as in the previous sections and we do not repeat all the computations 
but we only point out the main differences. It is worth to notice that this case is even easier to handle. 
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In particular, the irrelevant linear degree of freedom in the cost L(u) disappears: as explained below, the cost 
is now determined up to a multiplication by a positive scalar only and a simpler normalization is now given in 
Lemma 6.2. 

As for system (2.1), we consider the problems (Pz) and (44), but for the general system (6.1) now. 

Similarly to Section 3.1.2, we define the Hamiltonian function A(X, p, q,u) = pFo(q) + upF1(q) + AL(u), with 
p € R" and À < 0, in order to apply the PMP to Problem (Pr). 

Again, abnormal extremals have to be avoided, and we consider only normal extremals only, i.e., we assume 
that À = —1. 

All the information is still contained in the system: 


4 = Fo(q) + uF\(q) 
(6.2) 
L'(pF1(g)) + pFo(q) = 0. 


Using the u-parametrization of the experiments and taking into account the fact that p satisfies a linear ODE, 
the last equation of the previous system rewrites 


L(u) — uL'(u) = poV(u), (6.3) 
where po is the initial value of p(-), and V(u) is defined by 


poV(u) = p(t(u))Fo(q(é(u))). 


Most of the material and results developed in the previous sections remain unchanged with system (1.3) 
replaced by system (6.1). The main change is that the criterion L{u) is now reconstructed modulo multiplication 
by a scalar only. Indeed, the transformation L — aL(u) + bu does not leave invariant the set of extremal 
trajectories, only the transformation L — aL does. For this reason, Lemmas 3.14 and 3.15 modify as follows. 


Lemma 6.1. Let a be a positive number and let b be real. The two optimal control problems (Pr) and (Par) 
have the same set of extremal trajectories, while, in general, the two optimal control problems (Pr) and (Pr+vu) 
do not. 


Lemma 6.2. At uo, the cost function L to be reconstructed can be normalized so that L'(uo) = 1. 


Besides this, Theorem 4.1 is also true for system (6.1) and its proof generalizes straightforwardly, with some 
extra work. 
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Chapitre 8 


Conclusion 


Le travail effectué dans cette partie concerne la recherche d’une méthode anthropomorphique 
pour améliorer les trajectoires obtenues en utilisant les algorithmes de la partie I. 
En particulier, nous avons analysé le problème de contrôle optimal inverse associé aux drones. 


Le chapitre 6 introduit quelques méthodes bio-inspirées utilisées pour appréhender le comporte- 
ment d'animaux et définit le problème de contrôle optimal inverse visant à expliquer le comportement 
humain. 


Nous avons analysé le problème de contrôle optimal inverse pour des pilotes dans le chapitre 
7. Nous avons démontré que, pour un système général non-linéaire à un contrôle, il faut au moins 
trois expériences d’un même pilote pour reconstruire le coût minimisé. Dans le cas d’un drone de 
type HALE, le problème est légèrement différent et seulement deux expériences sont nécessaires à 
la reconstruction du coût minimisé. 

Ces coûts reconstruits sont utilisables dans une méthode de planification optimale pour construire 
des trajectoires proches de celles produites par des pilotes expérimentés. 
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CONCLUSION GÉNÉRALE 


Dans ce travail, nous avons étudié différentes méthodes de planification pour drones, de type 


HALE et MALE. 


La partie [ nous a d’abord permis de faire l’état de l’art sur les méthodes de résolution actuelles. 
Nous avons étudié dans le chapitre 2 les modèles cinématiques des véhicules considérés. À la fin de 
ce chapitre, nous avons présenté différentes méthodes de planification point-point, dans un milieu 
encombré ou non. 

De nouvelles méthodes de planification, permettant de répondre au problème de couplage vec- 
teur/capteur, ont été proposées dans le chapitre 3. En particulier, nous avons suggéré une méthode 
quadratique permettant d’asservir de façon pondérée la trajectoire du couple drone/caméra. 

Une de nos conclusions, dans le cas d’une caméra non contrainte, est que le problème de cou- 
plage vecteur/capteur se réduit au problème de planification du vecteur. Dans ce cas, une méthode 
indépendante de contrôle de la caméra a été proposée. 

Le problème de planification point-pattern a ensuite été traité. Deux méthodes ont été étudiées, 
l’une basée sur le principe de LaSalle et l’autre sur le PMP. 

Pour terminer, nous avons utilisé toutes les méthodes présentées pour proposer et appliquer un 
algorithme de suivi de cible mobile, en milieu contraint. 


Dans la seconde partie, nous avons étudié une méthode anthropomorphique pour améliorer les 
algorithmes de planification exposés dans la partie I. 

En particulier, nous nous sommes intéressés à un problème de contrôle optimal inverse et avons 

obtenu des résultats nouveaux, théoriques mais constructifs : 

e Au chapitre 6, nous avons introduit les méthodes bio-inspirées et le problème de contrôle 
optimal inverse. 

e Une méthode d’estimation du coût optimisé par les pilotes a été détaillée au chapitre 7. Nous 
avons prouvé que, pour un système général non-linéaire à un contrôle, il faut avoir au moins 
trois expériences d’un même pilote pour reconstruire le coût minimisé. Dans le cas d’un drone 
de type HALE, le problème est légèrement différent et seulement deux expériences sont néces- 
saires à la reconstruction du coût minimisé. 


Quelques perspectives pour continuer ce travail sont listées ci-dessous. 
Appliquer les méthodes de planification sur un drone 


Les méthodes étudiées au cours de ce travail ont été testées en simulation uniquement. Il serait 
intéressant de les tester dans un cadre plus réaliste. 


Réaliser la synthèse optimale en temps pour n’importe quel type de pattern 
La contribution que nous avons apportée concerne, entre autres, certaines méthodes de planifica- 


tion vers un cercle de rayon le rayon de courbure minimal du vecteur. Il serait intéressant d’étudier 
ces méthodes de planification pour tous types de patterns. 
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Étudier les méthodes de localisation de la cible 


Dans ce travail, la position de la cible est supposée connue à chaque instant. Cette hypothèse 
n’est pas toujours satisfaite. [Il faut donc mettre en place des algorithmes robustes de traitement 
d'images (ou toute autre méthode permettant de gérer des données capteurs) pour déterminer la 
position de la cible. 


Étudier le problème de contrôle optimal inverse associé aux drones MALE 


La partie IT présente des résultats concernant l’apprentissage à partir de données humaines. Nous 
avons considéré que les expériences faites vérifient un système de type Dubins. Cette hypothèse est 
valable uniquement si le drone étudié est de type HALE. Nous pourrions nous intéresser aux drones 
de type MALE. Pour cela, il faudrait étudier le problème de contrôle optimal inverse associé à un 
système différentiel non-linéaire à plusieurs contrôles. 


Étudier le problème de contrôle optimal inverse pour les drones HALE et MALE dans 
le cas où le pilote doit rejoindre un pattern 


Le problème de contrôle optimal inverse étudié correspond à des expériences de planification 
point-point. Or, nous avons vu que le problème de planification point-pattern est important. Une 
idée serait alors de considérer le problème inverse associé au problème point-pattern qui ajouterait 
une condition de transversalité sur la cible associée au problème de contrôle optimal. 
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Annexe À 


Article de conférence concernant le 
simulateur [6| : “Simulation of a UAV 
ground control station” 


L'article présenté a été accepté à la conférence MOSIM'12 : 9* International Conference of 
Modeling, Optimization and Simulation [6]. 
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ABSTRACT: In this article we present the development of a UAV ground control station simulator. We 


propose a module based description of the architecture of this simulator. 
a fixed wing aircraft. Finally we outline ideas for improved planification tasks. 
gures of the resulting station are displayed. 


through several diagrams, 


We recall the nonlinear model of 
The approach is made clear 


KEY WORDS: planification, simulation, UAV, optimal control, inverse optimal control problem 


1 INTRODUCTION 


Today’s interest in drone technology led to an in- 
creasing number of dedicated projects (Tisdale, Kim 
& Hedric 2009, Kim, Shim & Sastry 2002, Fabiani, 
Fuertes, Piquereau, Mampey & Teichteil-Künigsbuch 
2007, Ippolito, Yeh & Campbell 2009).  Those 
projects tackle problems such as the autonomy im- 
provement, the reduction of drone crashes due to poor 
availability of information, the organization of drone 
swarms or the calculation of secure and optimal flight 
paths. LSIS laboratory joins this research effort in the 
framework of project SHARE. Supported by a consor- 
tium of companies and research labs!, the research 
aims to develop à universal and interoperable ground 
control station for fixed and rotary wing UAVSs with 
a reduced number of operators. Specifically LSIS lab 
focuses on the connections between the UAV trajec- 
tory and its sensors. As a consequence, improved 
path planning algorithms that take into account pay- 
load requirements, optimal costs and obstacles (or no 
flight zones) avoidance are needed. 


In order to test and demonstrate our techniques à 
ground control station simulator is required. The 
use of engineering and simulation softwares allows to 
shorten development time and still maintain porta- 
bility of the code. Since the development is done 
in close contact with research partners, modular- 
ity is a key factor. Indeed, we want to be able 
to easily add, remove or modify parts of the sim- 


1Opéra Ergonomie, ONERA, Thales Alénia Space, Euro- 
copter, Adetel group. 


ulator. Another challenge is to determine the au- 
tomation level of the station. (Cummings, Platts & 
Sulmistras 2006, Sheridan, Verplank & Brooks 1978) 
propose insights on automation strategies that are 
useful to describe the simulator. On the Sheridan- 
Verplank scale, the station ranks no more than 3 and, 
according to (Cummings et al. 2006), it has a level 
of interoperability of 4 with respect to the STANAG 
4586 classification. This latter level can be described 
as follows: the ground station allows the control and 
monitoring of the UAV at the exception of launch and 
recovery situations. The simulation of such a device 
starts with the simulation of the aircraft’s trajectory. 
Then follows the emulation of all the data exchanged 
through the station between the operator and the 
UAV, and finally of the control algorithms. 


Classic tools for the implementation of the aircraft’s 
dynamics, the controller part and the man machine 
interface are Matlab and Matlab/Simulink. We pref- 
ered the use of S-functions to a systematic block based 
implementation in order to ease portability in case of 
discard of Matlab’s automatic code generation fea- 
tures. The simulation of the virtual environment 
and of the video flow is performed by à flight sim- 
ulator (Craighead, Murphy, Burke & Goldiez 2007). 
Flightgear which is open source and interfaces nicely 
with Matlab/Simulink is used (jun Yang, hui Qi & lin 
Shan 2009, Sorton & Hammaker 2005). 


The rest of this article is organised as follows. The 
structure of the simulator is presented in section 2. 
It is broken into several modules which functionali- 
ties are explained. The aircraft dynamics are recalled 
in section 3. Although project SHARE aims at a 
ground control station that addresses both fixed and 
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rotative wing aircraîts, we focus on fixed wing ones. 
Finally section 4 contains both low and high level de- 
scriptions of the control module (cf. figure 4). The 
reader is advised to reserve a special attention to sub- 
section 4.2 where important ideas regarding high level 
planification and learning algorithms are sketched. 


2 SIMULATOR 


A ground control station simulation instance has two 
main parts: the initialization phase and the scenario 
realization. The initialization phase allows to define 
parameters that are specific to a given simulation in- 
stance, this step is detailed in subsection 2.1. During 
simulation several different missions can be realized 
accordingly to the operator’s scenario. The architec- 
ture of this scenario realization part is detailed in sub- 
section 2.2. 


2.1 Simulator’s Initialization Phase 


We divide simulation parameters into two categories: 
setup parameters and core parameters. Core param- 
eters contribute to the simulation inner mechanisms: 
model parameters (subsections 3.3.1, 3.3.2 and 3.4) 
and control parameters (section 4). As such the user 
cannot manipulate them directly: they are either 
loaded through a setup assistant utility or managed 
through a modify/create assistant. 


The iînitialization procedure is schematized in figure 
1 below. The user selects a UAV type between fixed 
and rotary wings aircrafts which determines the equa- 
tions of motion that are solved. The model choice al- 
lows to specify which particular aircraft is simulated. 
The aerodynamic model (subsection 3.3.2), the corre- 
sponding aerodynamic parameters and the level one 
controller coefficients are loaded accordingly. Those 
can be further modified by the level one controls effi- 
ciency parameter (i.e. regular or degraded mode). 


Several payloads are embedded on the UAV. They are 
selected out of a list. In particular camera parameters 
are specified (subsection 3.4). Noise levels and bias 
of the sensors are also defined. High level planifica- 
tion algorithms constraints are specified in the form 
of minimum/maximum time to next waypoint, energy 
consumption related constraints, stealth mode, line of 
sight preservation or, more technically, weighting pa- 
rameters for the algorithms cost functions. 


Finally the set of maps of the area where the scenario 
takes place and the initial position of the aircraft? 
are selected. At the end of the initialization phase 
a Map synchronization utility ensures the connection 
between the 3D virtual environment and the 2D maps 
of the area of interest. 


2i.e. the origin of the world frame defined in subsection 3.2 
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Functional diagram of the initialization 


2.2 The Scenario 


In this subsection we propose a module based presen- 
tation of the simulator. This approach is schematized 
in figure 2. 
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Figure 2: Module based architecture diagram 


Let us start with the resolution of the UAV equations 
of motion. As said before they are completely deter- 
mined by the choîces made in the initialization phase. 
Input and output vectors are defined in subsection 3.1 
for a fixed wing aircraft. Let us focus on the output 
of this module, namely the state vector. It is used to 
generate the sensors outputs which are lumped into 
two categories: all sensors except the camera and the 
camera. The state vector is actually fully available 
from the several sensors à UAV carries. However we 
want to be able to simulate desynchronized sensors, 
noise levels and bias. The position which is available 
in the (x,y,z) system of coordinates is transformed 
into the WGS84 system for geo-localization purposes 
(Grewal, Weill & Andrews 2007). The UAV's atti- 
tude transforms from the unit quaternion representa- 
tion to the Roll-Pitch-Vaw angles one. The camera 
has a special treatment since first, a set of differential 
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equations has to be solved and second, the video flow 
is not generated in this module. We actually only 
take care of the line of sight orientation calculation 
and the evolution of extra camera parameters if they 
are used (subsection 3.4). 


The outputs of the two sensors related modules are 
sent to the geo-referencing module. This module en- 
sures the synchronization between the 2D map, the 
3D environment and the radar representation. The 
2D synchronization allows to display the trajectory of 
the aircraft on a map, to represent the video cone and 
the targets whose positions are known to the ground 
station. The radar gives information of the relative 
positions of possible targets. Finally the 3D synchro- 
nization is the computation of the actual video flow 
(ie. Flightgear). Those data are sent to the man- 
machine interface (MMI) module for display (figure 
3). This module is detailed further below in the 
present subsection. 


Al the modules described before feed the control 
module which in return sends its data to the equa- 
tions of motion and camera modules. It is the object 
of section 4. 


The exterior conditions module is the last part of this 
diagram, it encapsulates all the other ones since per- 
turbations can be introduced in all the other blocks. 
At this stage of the development we only use it to 
manage weather conditions, more particularly wind. 
À Von Karman model is used to generate the wind 
velocity vector required by the resolution of the equa- 
tions of motion. In our opinion, at least two extra 
modules shall be considered. We didn't add them 
in the diagram for clarity purposes and since they are 
at early development stages. Those are the extra pay- 
loads management (i.e. sensors specific to à mission, 
packages needed to be dropped.….) and the targets 
management unit. 


We end this subsection with a few words on the man- 
machine interface module. Figure 3 gives an account 
of how data flows from the operator to the UAV and 
back. Whereas the information sent back to the op- 
erator talk for themselves, it is not the case of all the 
features managed by the operator. UAV and payloads 
control modes are given in the next section. Simula- 
tion management features are start, pause, stop and 
visualization options. As for the rest, it is pretty ob- 
vious except for the mission type management part. 
This task is performed by the MMI modules which 
triggers the adequate module functionalities accord- 
ingly to the operator’s will. Let us focus on possible 
missions. The UAV has to perform several types of 
mission (NATO 2008) of which we implemented the 
following ones. 


Follow flight plan This is the basic mission. A 
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Figure 3: Man-machine interface data flow diagram 


flight plan is a series of waypoints the UAV has to 
visit. The flight plan is either defined before the 
scenario or made by the operator as the scenario 
undergoes. 


Patterns This mission consists in following operator 
parameterized patterns: circle, racetrack, and 
figure eight (table 1) in order to detect any ac- 
tivity in a specific area and/or collect informa- 
tion. The trajectory that leads to à given pat- 
tern is the subject of more or less sophisticated 
algorithms (section 4). Indeed the UAV has to 
be tangent when it enters the pattern. The ap- 
proach trajectory to the pattern is constrained 
by the geographical configuration, UAV physical 
characteristics, UAV inital orientation and the 
traveling direction of the pattern. 


Tracking The objective is to follow a friend or en- 
emy point of interest. In the case of an ally, an 
expanded surveillance zone is considered. Mul- 
tiple situations can be encountered: the target 
is in à no-fly zone, target occlusion occurs, the 
target is moving, camera commands are either 
automated or manually activated by the oper- 
ator. In the case of an enemy target, it is of 
utmost importance to be able to maintain it in 
the camera field of view. 


Replanification Replanification happens when the 
UAV needs a waypoints list that is better gener- 
ated with appropriate algorithms — i.e. when a 
mission has to be changed because of unexpected 
events. No-fly zone avoidance, search of an ad- 
missible path and optimization of an optimality 
criteria are the three aspects that drive this mis- 
sion type (subsection 4.2). 
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pattern parameters list 


common to all center point 


rotation direction 


circle radius 
racetrack / small radius 
hippodrome large radius 
orientation 
figure eight pseudo radius 
orientation 


Table 1: Pattern parameterization 


3 MODELING 


3.1 Notations and Definitions 


u() vector v projected in frame (f) 
Id (n x n) dimensional identity matrix 
X() aircraft’s position vector 
vU) aircraft’s speed vector 
CP ) aerodynamic speed vector 
q unit attitude quaternion 


matrix associated to vu + Q x v 


aircraft’s (resp. camera’s) roll angle 


O(%) (resp. 4e) aircraft’s 


y(w) (resp. we) aircraft’s (resp. cameras) yaw angle 
uw (A) aircraft’s angular velocity vector 
(Le, Yes Ze)(®) camera’s position 

a aerodynamic angle of attack 

B aerodynamic sideslip angle 

F0) aerodynamic force 

FE propulsive force 

M, aerodynamic moment at point g 
M, propulsive moment at point g 

(e162) gravity vector (ie. G(%) = (0,0,g)) 
I aircraft’s inertia matrix 

on roll control 

dn pitch control 

Ôn yaw control 

Ôx thrust control 

T' overall thrust 

OPA pe) set point 

Ô9 gt) set point 

OM Lo) set point 

u(t) (resp. uc(t))  UAV (resp. camera) control vectors 
p atmospheric density 

Po atmospheric density at sea level 

S aircraft’s surface of reference 

l aircraft’s length of reference 

m mass of the aircraft 


The controls are 


2012 - 


(resp. camera’s) pitch angle 


Bordeaux - France 


and the state vector is 


T 
qg = id . 
0) 0 y 


Je 


The matrices LG) . and Î are defined as 


: 0 —qg3 q 
[(G) | = q3 0 -q 
a  gq 0 
— 92 1 


and 1 = IT = (is € {1,2,3}°) 


3.2 Hypotheses and Frames 


The model for aircraft dynamics proposed here is in- 
spired by (Boïffier 1998, Wanner 1984, Junkins & 
Schaub 2001). The modeling hypotheses are° 


1. earth is flat and motionless with respect to the 


simulation duration, 


2. gravity is constant (i.e. the center of mass is the 
same as the center of gravity), 


. atmospheric pressure and temperature are con- 
stant (i.e. only density variation with respect to 


height is considered), 


. there is a uniform wind velocity field, 


Qt 


. the aircraft is a 6 DoF solid which fuselage has a 
symmetry plane, 


. the propulsion system is on the fuselage axis, 


. the mass and the inertia matrix of the aircraft 
are constant parameters. 


The first hypothesis makes the world frame inertial, 
thus we only need three frames to derive the equa- 
tions. 


World frame (w) This frame is attached to a refer- 
ence point ©, the Ox, Oy, Oz axes are oriented 
in the north-east-down directions. Although the 
aircraft’s position expressed through the equa- 
tions of motion is given in the (x, y, z) system, we 
need to express it in the WGS84 system in order 
to perform geo-referencing (Grewal et al. 2007). 


3See (Boiffier 1998) for a comprehensive treatment of the 
hypotheses. 
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Body frame (b) This frame is attached to the air- 
craft’s center of mass (G). The Gx) axis goes 
along the fuselage axis. The G2() axis is taken in 
the plane of symmetry of the aircraft and points 
downward, and Oytt) = Oz) x Oxtt) 

The composition of a translation and a rotation 
transforms frame (w) into frame (b). The ro- 
tation defines the attitude of the aircraft. We 
use the Euler angles Roll-Pitch-Yaw parameter- 
ization for its physical meaning, and the unit 
quaternion one to solve differential equations 


(Junkins & Schaub 2001). 


Aerodynamic frame (a) This frame shares his ori- 
gin with frame (b). The Gx(%) axis is carried by 
the aerodynamic velocity vector. The rotation 
that transforms Gx() into Gx(®) is parameter- 
ized by the angle of attack and the sideslip an- 
gle. The other axes are obtained through this 
rotation (subsection 3.3.2). 


3.3 Aircraft Dynamics 


3.3.1 Equations of Motion 


The dynamic resultant theorem and the dynamic mo- 
ment theorem write as 


wE) x mvE) 
(1) 
160) = MP, + MP, - 00 x Iu®) (2) 


mEV® = FO + FO + MED 


e,g 
where 
() (b) 
po) = PSV e FO 2e 
2 se J ot 
—CY Forz 
(3) 
() 
C: Ma 
SIV2 ! 
M) =, À ; ae cs M) _ M 
Ca M 


The rotation matrix from the world frame to the body 
frame expressed in terms of q, and the dynamics of q 
are (Junkins & Schaub 2001) 


R(g) = (aÿ — 44") Id3+2(44-&[g"]) (5) 


and 


1 0 710) 
= | 6 [CON q (6) 
Therefore we have 


mG) = R(q) (0,0,mg)" (7) 
X(w) = le) = RT(q)V®) (8) 


The equations of motion of the aircraft with state 
(x), 0), q,w)) are equations (8), (1), (6) and 
(2). In order to actually use them, we still need to 
define à few elements. 


e The atmospheric density model is 


—1.1210 4h 


—4 
p = pot = pper 1210 "2 


e From hypothesis (5), elements Z12 and Z3 of 
the inertial matrix are null (the symmetry of the 
inertial matrix implies 121 = 132 = 0). 


e From hypothesis (6), the propulsion force vec- 

tor components F9), and FE, are null. The 
propulsion moment vanishes — see (4.3.4) in 
(Boïffier 1998), with @n — Bm = 0. 
An engine efficiency model can be used for the 
modulus of the thrust (1e. F®} in our case). 
For example, in (Boiffier 1998) T = k»PVA ôx 
is proposed. We simply used T = k,,0%, with 
O0 < km < 1. 


e The wind is characterized by its velocity vec- 
tor pe). The aerodynamic speed vector is then 
defined as v@) = v:®) — V4), and V2, = 


(VE) ve 


3.3.2  Aerodynamic Model 


We call aerodynamic model the equations used to 
fully express aerodynamic forces and moments (equa- 
tion 4). Force coefficients are first defined in frame 
(a) (.e. from wind tunnel experiments), then rotated 
into frame (b). As such the angle of attack and the 
sideslip angle are needed. Starting from the aerody- 
namic speed vector we derive equations 


ve, = Vie cos(æ) cos(B) 


VO) = Vs sin(8) 


ae,y 


VO) = VV, sin(a) cos(B) 


ae,Zz 


and we write 


26® sin(@) 


C1 c®) cos(æ) cos(B) + cf sin(B) cos(a) 
cc = cf cos(B) — c® sin(B) 
el c® sin(a) cos(B) + cf sin(B)sin(a) 


+c® cos(a) 


There are many possible models for cé), ch), GER 


CYy, Cm and Cn (Wanner 1984, Schmidt 1998). We 
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cite the simple model 


Cire Core 
ae = Cysl + Cu, ôn 
= Corte 
G = Capsin(s) +1 (Cie, Se + Crus æ) 
+ C6, di + Ci,5, ê 
Cm = Cas FOnée dns + Crus ge 
Cn — CB Le Cnôn Ôn + Cds Ve 


3.4 Camera Dynamics 


This model tracks the orientation of the line of sight 
of the camera. This latter is assumed to be solidly 
attached to the center of mass of the aircraît: 


rt = Xx{v) 
y = X{) 
ZX ce 


Following camera datasheet specifications (wescam 
2011) the dynamics can be approximated as a first 
order response to an attitude set point. This atti- 
tude set point is given in the (b) frame (i.e. pe) = 
(68) + 60) /T). The time constant 7 defines the 
speed of the camera steering system. Constraints on 
the coverage can be expressed in terms of constraints 
on the inputs rather than on the output (a significant 
asset both for implementation and optimal control 
purposes). 


60 = nf) — 1 (GE — gtu)) + La 
jp) = 10 om) +15 (9) 
je = 1 (9 yo) 1, 


Extra parameters that can used to model the camera 
system are described in the following table (wescam 
2011). 


parameter indicative value 
optical eye field of view 36 to 1 (°) 
IR fields of view 30, 7 and L.8 (°) 
frequency 24 (FPS) 
turret coverage 

azimuth O0 to 360 (°) 
elevation 90 to -120(°) 
roll na 


turret steering speed 0 to 60 (°/sec) 
turret stabilization quality  abt 3.107% (°) 
autofocus time constant 0 

(ie. instantaneous) 


4 CONTROL MODULE 


As explained in (Cummings et al. 2006) both the level 
decomposition and strategies of automation are cru- 
cial elements in order to achieve efficient accurate and 


safe UAV operations. First we give the global picture 
of the control module in subsection 4.1, and then fo- 
cus on two high-level aspects in subsection 4.2. Those 
two important topics are planification and learning 
techniques. 


4.1 General Considerations 


Three levels are considered for this module (figure 4). 
Level zero corresponds to a direct access to the con- 
trol surfaces of the UAV. Although unrealistic we kept 
the possibility to directly manipulate control surfaces 
as far as the simulator only is concerned. Level one 
takes an input signal of the form heading, speed, alti- 
tude and outputs à (ô, dm, 0,04) vector (section 3). 
This controller consists in a stability augmentation 
system (yaw, pitch, and phugoid dampers) in series 
with à basic PID like autopilot system (heading, air- 
speed and altitude holders). Note that the complex- 
ity of this controller grows with the complexity of the 
aerodynamie model (recall that the controller param- 
eters are loaded accordingly to the UAV model, cf. 
section 2). Finally level two encapsulates all the al- 
gorithms that generate à heading, speed, altitude set 
point. 


jme mm mme & 


Ne 
\ 
A 
REPLANFICATION 
ALGORITHM à 


LEVEL 2.3 


LEVEL 2.2 


GALCULATION 
LEVEL 2.1 \ OF NEXT 
HEADING / SPEED / ALTITUDE 


HEADING / SPEED / ALTITUDE 


PNA CONTROLLER 


ÿ 


UAV 
LME CONTROLS 


Figure 4: Level based decomposition of the controller 
module 


The operator has the choice between several different 
piloting modes (figure 3) for both the UAV and the 
camera turret. 


Manual mode Through a joystick, the heading, 
speed, altitude setpoint is directly controlled. As 
mentioned before, a very special option, only 
available in simulation, grants access to the con- 
trol surfaces of the UAV. 


Semi-manual mode At least one of the three com- 
ponents of the heading, speed, altitude setpoint 
is managed by the system — e.g. altitude is kept 
fixed. 
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Automatic mode The level one set point is defined 
by a level two algorithm which depends on the 
active mission. 


The control modes of the camera are simpler to han- 
dle as only manual and automatic modes are pro- 
posed. Knowing the lastly visited, the targeted and 
the next on the list waypoints, the calculation of 
the new heading, speed, altitude set point is easily 
obtained from geometry. Note that a waypoint is 
defined by its WGS84 coordinates and an optional 
timestamp. When minimum/maximum time between 
waypoints constraints are activated, timestamps are 
used to adjust the set speed. 


The next three parts to explain are the pattern, re- 
planification and pursuit algorithms. As said before à 
pattern is defined by a set of parameters (table 1) as 
such a waypoints list describes it. What remains to be 
done is to propose à trajectory that allows the UAV 
to reach the pattern with the correct configuration 
(ie. tangentially), starting from any initial point and 
any initial orientation. This is one of the jobs done 
by the replanification algorithm. The second one is 
to generate optimal trajectory to travel through an 
area without having a pre-established flight plan. A 
strategy in order to use this algorithm is to solve a re- 
planification problem, then propose the solution tra- 
jectory as a list of waypoints and finally refresh the 
list through a new calculation each time à waypoint is 
reached. The pursuit algorithm output can be seen as 
a heading, altitude, speed vector that is recomputed 
whenever necessary. 


À few insights on replanification algorithms are pro- 
posed in the next subsection. We even make one step 
farther by presenting experiment based learning of 
optimal cost functions. 


4.2 Planification and Learning 


4.2.1 Introduction 


In the context of the SHARE project we have also to 
achieve two extra tasks: 


1. Fill in the module called planification algorithms 
in figure 4. There is a lot of bibliography on this 
topic (Betts 2001, Bullo & Lewis 2004, Fliess, 
Lévine, Martin & Rouchon 1995, Kim et al. 2002, 
Laumond 1998, LaValle 2006, Park, Deyst & 
How 2004, Van Nieuwstadt & Murray 1998) , 
non exhaustively. These classical algorithms are 
based upon very different approaches such as ge- 
ometric control, optimal control, flatness. In fact 
we are also developing our own methods that we 
present briefly below (subsection 4.2.2). 


2. An idea to develop planification/replanification 


methods is to learn from the behavior of experi- 
mented pilots. We do this here for HALE (High 
Altitude Long Endurance) drones. In subsection 
4.2.3 we present briefly our ideas, that are in- 
spired from the beautiful work of Jean and al. 
in the papers (Chitour, Jean & Mason to be 
published, Chittaro, Jean & Mason to be pub- 
lished). In these papers, they attack the prob- 
lem of identification of the cost minimized in hu- 
man locomotion. Our problem is very similar 
since HALE drones behave kinematically more or 
less as a human being moving on à plane (con- 
stant altitude, constant speed). There is a lot 
of other methods dedicated to this human loco- 
motion problem, see for instance (Li, Todorov & 
Liu 2011, Berret, Darlot, Jean, Pozzo, Papaxan- 
this & Gauthier 2008). 


4.2.2 Planification/Replanification 


In general the mission of a drone is planed in advance 
and specified by à certain number of checkpoints and 
a certain number of patterns (line, circle, figure eight 
and hippodrome). There is the need of on-line replan- 
ification methods when the mission is interrupted and 
the drone has to: 


e Join a fixed target and turn around following à 
certain pattern, 


e Join à moving target and follow it. This is not 
an obvious problem when the minimum speed of 
the drone is higher than the target speed. 


Regarding the case (not yet treated above) of HALE 
rotorcraft-based drones, they behave kinematically 
more or less as the classical “simple car” model 
(Laumond 1998). Hence one can formulate an opti- 
mal control problem which looks like a left-invariant 
subriemannian problem over the group of motions of 
the plane. On this topic there is the beautiful com- 
plete mathematical work of Yuri Sachkov (Moiseev 
& Sachkov 2010), Sachkov 2010, Sachkov 2011) that 
does the job, if no obstacles are taken into account. 
If obstacles occur, the Sachkov method can be cou- 
pled to a method for finding first admissible trajec- 
tories avoiding obstacles and approximating them by 
pieces from the Sachkov synthesis. Several methods 
are available to find such admissible trajectories see 
(LaValle 2006). 


The case of fixed wings drones is different. If we 
consider à simple model of Dubins type (Laumond 
1998), the problem can be stated as à repetition of 
minimum-time problems for the Dubins car, where 
the final target is a (non-oriented or oriented) pat- 
tern. This optimal control problem can be solved eas- 
ily, but it leads to à non-smooth and eventually non- 
continuous optimal synthesis. In fact, a nice smooth 
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optimal synthesis can be found, leading to trajecto- 
ries shown in figure 5. Details on this optimal syn- 
thesis and proof of the stability and convergence will 
be given in a forthcoming paper. 


Figure 5: Trajectory resulting of smooth optimal syn- 
thesis starting at the point (20;0) with a direction of 
0 radian and arriving tangent to a pattern, here the 
circle centering at (0;0) 


4.2.3 Learning from Experimented Pilots 


When the planning problems are specified in terms of 
optimal control problems, the choice of the cost to be 
minimized can be made from the experience of exper- 
imented pilots. For our fixed wings HALE drones, as 
we said, the basic kinematic model is the classical Du- 
bins one. It has been shown by Jean and al. (Chitour 
et al. to be published, Chittaro et al. to be published) 
that the cost to be minimized for human locomotion 
is an integral length-curvature compromise. We com- 
pleted their work by developing à program allowing 
to identify this cost. Moreover, we proved the follow- 
ing general theoretical result: for a generic model of 
motion, the integral cost can be recovered from two 
experiments, provided that these experiments visit at 
least twice the same value of the control. These de- 
velopments will be the purpose of another paper. On 
the figure 6, we show the reconstruction of the cost 
as à function of length and curvature from two exper- 
iments. 


5 CONCLUSION 


In this article we presented a ground control station 
simulator for fixed wing aircrafts that can easily be 
extended to rotative wing ones. This description has 
been done in the form of a module based decomposi- 
tion. The control module was given special attention 
as it deals with several issues important in order to 
improve the UAV autonomy. We sketched ideas to 
tackle the planification/re-planication task which can 
be used to deal with pattern based navigation and 
tracking problems. We also presented an approach to 
improve the cost functions that are used in optimal 
control synthesis based on trajectories obtained from 
piloted aircraîts. 
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Figure 6: Reconstruction of the cost as the function 
of the curvature and the theorical cost 


To conclude this work we display in figure 8 a view of 
the ground control station simulator. Snapshots of a 
monitoring mission of a famous building while flying 
in circles are shown in figure 7 (white squares have 
been added to ease visualization). 


Figure 7: Keeping a close watch on the tower 
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Figure 8: Large view of the simulator 
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Annexe B 


Article de conférence concernant 
l’application des méthodes de 
planification sur le simulateur [4] : “Path 


planning and ground control station 
simulator for UAV” 


L'article présenté a été accepté à la conférence ZEEE Aerospace Conference ||. 
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Abstract—In this paper we present a Universal and Interoper- 
able Ground Control Station (UIGCS) simulator for f xed and 
rotary wing Unmanned Aerial Vehicles (UAVSs), and all types of 
payloads. One of the major constraints is to operate and man- 
age multiple legacy and future UAVS, taking into account the 
compliance with NATO Combined/Joint Services Operational 
Environment (STANAG 4586). Another purpose of the station is 
to assign the UAV a certain degree of autonomy, via autonomous 
planif cation/replanif cation strategies. The paper is organized 
as follows. 


In Section 2, we describe the non-linear models of the f xed and 
rotary wing UAVSs that we use in the simulator. 


In Section 3, we describe the simulator architecture, which is 
based upon interacting modules programmed independently. 
This simulator is linked with an open source fight simulator, to 
simulate the video f ow and the moving target in 3D. To conclude 
this part, we tackle brief y the problem of the Matlab/Simulink 
software connection (used to model the UAV”s dynamic) with the 
simulation of the virtual environment. 


Section 5 deals with the control module of a fight path of the 
UAW. The control system is divided into four distinct hierarchical 
layers: fight path, navigation controller, autopilot and fight 
control surfaces controller. 


In the Section 6, we focus on the trajectory planif ca- 
tion/replanif cation question for f xed wing UAV. Indeed, one of 
the goals of this work is to increase the autonomy of the UAV. We 
propose two types of algorithms, based upon 1) the methods of 
the tangent and 2) an original Lyapunov-type method. These 
algorithms allow either to join a f xed pattern or to track a 
moving target. 


Finally, Section 7 presents simulation results obtained on our 
simulator, concerning a rather complicated scenario of mission. 
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1. INTRODUCTION 


Recently large technological advances have been realized in 
the felds of UAVS in order to reduce crashes and collisions 
[11, [2], [31]. Human error remains the main cause, as reported 
by the NATO Joint Capability Group for UAVSs. For this 
reason it is important to improve on the autonomy of UAVSs 
via the coupling between the carrier (UAV) and the payload 
(camera for instance). 

To meet these requirements, a strong consortium formed by 
several companies and research labs! has been created. This 
consortium has decided to join the R&D efforts of each 
partner within a project called SHARE. 

SHARE project consists of developing a universal and in- 
teroperable ground control station for f xed and rotary wing 
UAVSs and all types of payloads, with a reduced number of 
operators at the level of the base station. 

Several technical requirements must be treated to achieve this 
project: 


1. Increase mission effectiveness by a) controlling the trajec- 
tory ofthe UAV and b) merging several work stations. 

2. Improve on data transmission by securing data links and 
reducing compression and decompression time of images. 

3. Increase the universality of the ground station to support 
any type of UAV and payload. 


Specif cally LSIS lab focuses on developing a simulator of 
a ground control station and developing algorithms for path 
planning. Those algorithms must take into account payload 
requirements, optimal costs and obstacles (or no fight zones). 
The use of fight simulation tools for complex aerospace sys- 
tems to reduce timetable, risk and required amount of fight 
testing is a well-recognized advantage of these approaches. 
For these reasons we have developed a ground control station 
simulator which will serve several purposes: 


. Simulate the behavior of the carrier, its onboard payloads 
and information transmitted to the control station. 

- Provide support for the simulation of different missions 
(tracking, monitoring, etc.). 

. Serve as a development platform to test several path plan- 
ning algorithms. 

. Serve as a demonstration platform. 


The simulator is developed under the Matlab/Simulink en- 
vironment. It simulates the non-linear dynamic equations 


Thales Alénia Space, Opéra Ergonomie, Eurocopter, ONERA, Adetel 
Group 


of UAV, payload dynamics, control module, MMI (Man 
Machine Interface), exterior conditions, etc. 


In Section 3 we present the simulator architecture based 
upon interacting modules. Graphical representations used in 
the simulator are discussed. In 3D, to represent the video 
fow and the moving target we use an open-source fight 
simulator: Flightgear. Otherwise, we have developed our 
own 2D simulation environment to visualize UAV”s trajec- 
tories. À synchronization between 2D and 3D environments 
1s established. 


In Section 5, we deal with the control module and its different 
layers. The control module calculates commands to send 
to the control surfaces of the carrier and the camera. The 
automation strategies of the control system allow achieving 
accurate and safe UAV operations. Four levels are considered 
for this module which offers ability to choose the control 
mode and the type of mission (collecting information, moni- 
toring, tracking, etc.). 


In Section 6 we present our solutions to the problem of 
trajectory planif cation/replanif cation for fxed wing UAV. 
We propose original algorithms based upon the tangent and 
Lyapunov vector f eld methods in order either to join a pattern 
or to track a moving target. 

Considering a fxed target, missions consist of joining and 
following operator’s parameterized patterns: circle, racetrack 
or leminscate in order to detect any activity in a specif c area 
and/or to collect information. 


Section 7 is devoted to simulation results for several scenarios 
in order to test and validate our simulator and replanif cation 
algorithms. 


In section 9, with our conclusions, we outline ideas for 
position estimation based upon knowledge of the target and 
environment information. Another idea for future is to de- 
velop a decision support system in order to aid the operator 
in a replanif cation context. 


2. UAV MODEL 


This section deals with the modeling of UAV”s fight dynam- 
ics mechanical equations. The model is the heart of the sim- 
ulation process; it describes the dynamic and aerodynamic 
behaviors of the carrier. 


For the sake of generality we try to keep generic equations. 
They correspond to two types of carrier: fxed and rotary 
wing. The simulation problems are similar in both cases how- 
ever the path planning problems are much more complicated 
in the fxed wing case. It is the reason why we focus on this 
case here. 

To simulate a specif c carrier it is enough to specify its aero- 
dynamic data in a way that will be described later (Section 
3). 
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Notations and Def nitions 


vu) 
XO) = (x,y,2) 
(D, 8, p)() 


do, 41, q2, 3 
y (4) = (u,v, w)() 


Mérep 

(ôl, 0m, ôn) 
Cx, CY, CZ 
C, Cm: Cn 


Xmrs Ymrs Zmr 
X fus) Yfus; Z'fus 
X4r, Ver, Tir 

T 

Lonr Mnr Nonr 
Lir, Nir 


Oco1,0 
(dtat ; Ôlon: Ôcot Ôped) 


Les Ye. ze) 


( 
(De, 6, Ye) 
(de, 06e, ôbe) 


: vector v projected in frame (f) 
: aircraft’s position vector 
: Euler angles of aircraft 


(roll, pitch, yaw) 


: unit attitude quaternion 

: aircraft’s speed vector 

: aerodynamic speed vector 

: aircraft’s angular velocity vector 

: aerodynamic force 

: propulsive force 

: aerodynamic moment 

: propulsive moment 

: aircraft’s commandés (roll, pitch, yaw) 
: aerodynamic force coeff cients 

: aerodynamic moment coeff cients 

: mass of the aircraft 

: gravity vector (i.e. gl“) = (0,0, g.)) 
: aircraft’s inertia matrix 

: atmospheric density 

: aircraft’s surface of reference 

: aircraft’s length of reference 

: aerodynamic angle of attack 

: aerodynamic sideslip angle 

: force applied to the helicopter 


: moment applied to the helicopter 
: aerodynamic forces generated 


by main rotor 


: aerodynamic forces generated 


by fuselage 


: tail rotor forces 
: main rotor thrust 
: aerodynamic moments generated 


by main rotor 


: aerodynamic moments generated 


by tail rotor 


: main rotor induced velocity 

: main rotor rotating speed 

: main blade radius 

: lift curve slope of the main rotor 


blade 


: main rotor blade number 
: main rotor blade chord length 
: longitudinal and lateral 


tip-path-plane (TPP) f apping angle 


: ratio of main rotor blade collective 


pitch to üco 


: main rotor spring constant 
: main rotor hub’s vertical position 


above the center of mass 


: prof 1 power of main rotor 

: induced power of main rotor 

: parasite power of main rotor 

: climbing power of main rotor 
: trim offset of the main blade’s 


collective pitch angle 


: helicopter’s commands 


(aileron, elevator, collective, rudder) 


: camera’s position 
: angles of camera (roll, elevation, azimuth) 
: camera”’s commands 


UAV Model for f xed wing 


Mainly, we use the dynamical equations from [4], [5], [6]. 
Basically they are just the Euler equations of the solid. 
Modeling assumptions are the following: 


1. The aircraft is a standard 6 DoF solid, symmetric w.r.t. its 
vertical middle plane 

2. The mass and the inertia matrix of the aircraft are constant 
3. The propulsion system lies along the fuselage axis 

4. Atmospheric pressure is constant and effects of external 
temperature are neglected 

5. We work in a fxed coordinate frame in which the move- 
ment and the roundness of the earth are ignored. This frame 
1s called the “world frame”. 


In fact, to describe the motion we need to consider three 
frames: 


World frame (w) This frame is attached to a reference point 
O. The Or, Oy, Oz axes are oriented in the north-east-down 
directions 

Body frame (b) This frame is attached to the aircraft’s cen- 
ter of mass G. The Gx(Ÿ) axis goes along the fuselage axis. 
The Gz() axis is taken in the plane of symmetry of the 
aircraft and points downward. Thus Oylt) = Oz) x Or) 

Aerodynamic frame (a) This frame shares its origin with 
frame (b). The Gx!{9) axis is carried by the velocity vector. 
The rotation that transforms Gx(*) into Gx(®) depends on the 
angle of attack and the sideslip angle. For precise def nition 
of the other aerodynamic coordinates see [4]. 


To avoid the singularities of Euler angles, we prefer to use 
unit quaternions in the computations. However Euler angles 
Roll-Pitch-Yaw parameterization is used to deliver informa- 
tion to the user [5]. 


We just show here the quaternionic rotation from frame (b) to 
(w): 


2(q193 + gog2) 
2(q2q3 — qoqi) 
2(g5 + gs) — 1 


R = | 2{(qiq2 + gog3) 2(Q8 + d) — 


2(q6 + qi) — 1 2(q1q2 — 403) 
2(qg1q3 — qog2) 2(q2q3 + qoqi) 


Finally, the state vector of our system has dimension 12: po- 
sition (3), velocity (3) , quaternions (3) and angular velocity 
(3). Due to the quaternionic representation we integrate the 
13-dimensional differential system (3), (4), (5) and (6). 


The six-degree-of-freedom rigid-body dynamics can be ex- 


pressed by the fundamental equations, for translation and 
rotation: 


mVO +00 x mV 0 = FD, + F0, +mg0 (D 


160) +00 x Iu0 2 MG, +MP, © 
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where 
FDo= | SPSVaCy 


| 2(q193 — Gog2)Mmgz 
mg) = 2(q2q3 + gogi)Mmgz 
(2 


Fr] Fe 
ne = |F,] = | 0 | (by applying assumption 4) 
F,] (ù 
| sci] 
9 aC1 
1 
Mièro = | 5PSVaCm 
- SC 
2P an 
b) 
M, = 0 
A —F -E 
I=| -F B -D 
-E D C 


With À,B,C,D,E real constants and F = D = 0 by 
applying assumptions 1 and 2. 


By applying the equation of forces (1) we deduce the accel- 
eration vector in the frame (b) 


; p Fr 
ù = ——SVÈCx + — — wow + wav + 2g,(q193 — Q0q2) 
2m m 


F 
L — Liu + wau + 2g,(g2q3 — Goq1) 


ù = L-SV2Cy + 
2m 


m 


: P 2 z 2 2 
= ——SVSCz + — + gz(2 — 1 
de Se SVa Oz + + — univ +wau + g:(2(q0 es ) 
According to the equation of moments (2) we obtain the 
angular acceleration vector in the frame (b): 


(AC — E?jin = Psv2(GiC POLE) 


+ (BC = C? EL E? wow — E(B — À — C'jwiwa 


l 
Bi = ESVÈ On OA Grue 


(AC — E?jùs = PSVÈ(CE 
+ (E? DE A(B = A))wiwo 


(4) 
The basic kinematic equation of the solid provides the motion 
of the quaternions: 


d2 —2| uw —ws ù w q2 6) 
d3 w3 Wa —Ww] 0 q3 
The aircraft position vector X (w), in the inertial frame, is 


given by 
X(&) = rRy) (6) 


Aerodynamic model 


The aerodynamic model consists of the data necessary to 
compute the external forces and moments. We have chosen 
to represent the aerodynamic data under the guise of lookup 
tables. The advantage of this choice is the representation of 
non-linearities. Other standard methods consist of deriving 
linear relation. In [7] the authors describe two methods and 
present the advantages of using lookup tables. It is possible to 
combine several approaches. In our simulator lookup tables 
are read and stored during the initialization step. 


More precisely, we consider a set of six lookup tables 
which correspond to the three aerodynamic force coeff cients 
Cx,Cy,Cz and the three aerodynamic moment coeff cients 
Ci, Cm; Cn. These coeff cients depend on angle of attack, 
sideslip angle and f ap defection. Equations of force coeff - 
cients are: 


Cx = Cro + HO 


Cy = CysB + Cyp V, EE 72 + Cysn ON 
wol 
Cz = Ca + Cia a + Crsm ÔM 


Equations of moment coeff cients are: 


l l 
Gi = Ci,8 + Ci, D + OT + Ci on + Ciuôl 


Ci = Cet On “ut La 
Ve 
l l 
Cn = Cons B + On . + Cnsn On + Cr ôl 


Coeff cients Cr, Cygs Co Clg Cmar Cng are also tabulated 
with respect to values of the f ap defections. 

Standard interpolation methods are used to extract the data 
from these tables. 


UAV Model for rotary wing (sketch) 


Our model for a rotorcraft (single main rotor with tail rotor) 
is detailed in [8], [9]. 

Again from the fundamental equations (1) and (2) we obtain 
a similar 12-dimensional differential system. Assumptions 
made in section 2 still hold true. 


The equation of forces and moments which are applied to the 
carrier are respectively given by: 


Fr Xi + X fus 
F, —= Fy LL Von L Vus L ir 
F2 Znr + Zfus 
Mix Lnr ca Lir 
M, = | My| = Mynr 
M Nr + Nir 


To compute the main rotor forces and moments, we frst 
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camera) 


Camera 


Targets 
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Figure 1. Module based architecture diagram 
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calculate the main rotor thrust T and the induced velocity v;: 
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Next, we compute the forces and moments generated by the 
main rotor as follows: 


Xmr = —Tsin(a;s) 

Ynr = Tsin(b:) 

Zmr = —T cos(a,) cos(b;) 

Lnr = (Kg + THyr) Sin(bs) 

Mr = (K8 + THmnr)sin(as) 

Nr = —(Ppr + Pi + Poa + Pe)/wmr 


The f nal differential system is given by (3), (4), (5), (6) with 
the previous values of forces and moments. 


3. SIMULATOR ARCHITECTURE 


The goal is to defne a completely modular architecture al- 
lowing operators to easily add, remove or modify some parts 
of the simulator while maintaining the real characteristics 
of a ground control station. The user can beneft from the 
help of assistants modules to add or modify the aerodynamic 
coeff cients and ”’payload” models. 


This system must be able to simulate all the operational 
environment of the station, in order to test and verify the 
development of this one. It should allow to validate our 
autonomous path planning algorithms. Furthermore, the 
simulator can be used to train operators or as a portable 
demonstrator. The architecture must take into account all 
these aspects. 


In this subsection, we brief y describe the different modules 
of the simulator. A schematic of the interactions between 
each module is shown in Fig.1. 


The module carrier ensures the resolution of the UAV'’s 
equations of motion described in section 2. The outputs of 


the controller and the exterior conditions are the inputs of 
this module. The output consists of the state vector, as a 
continuous function of time. 


The module Camera deals with the camera’s equations of 
motion. The model tracks the orientation of the camera’s line 
of sight. The camera is assumed to be rigidly attached to the 
center of mass of the aircraft. 

For the camera, we use the following simple f rst order model: 


Te = av) 
Ve — y) 
ze = 2%) 


: 1 1 

QE) = un — = (al — gr) + 2og, 
T x 

| 1 1 

LOEPRES (es _ ni) + 268 
T T 

- 1 : 1 

DEN = ua — 2 (ul y) + Loue 
HA T 


The variable 7 def nes the time constant of the camera steer- 
ing system. The input of this module is the state vector of 
the UAV, the controls, the target data and exterior conditions. 
The output is the state vector of camera. 


The Control module ensures the control of the UAV and the 
camera. Ît has several modes and levels of control. We detail 
this module in section 5. The input of this module is the sate 
vector of the UAV and the camera, the exterior conditions, 
the geo-referencing data and the set points given by the user 
via the MMI. The output provides the controls of carrier and 
camera. 


The Exterior conditions module manages weather conditions, 
more particularly wind. We use a Von Karman model [10] to 
generate the wind velocity vector required by the resolution 
of the equations of motion. Indeed, the wind speed V,, 
allows us to calculate aerodynamic speed vector V, and de- 
duce thereafter aerodynamic angles a and 5 from geometric 


relationships between the speed vectors V:9) and V0). 


MMI (Man Machine Interface) manages the fow of infor- 
mation between the operator and the UAV. This unit allows 
the user to control the carrier and the payload, and to def ne 
the mission to achieve (monitoring, collecting information, 
tracking, etc.). Also, it displays the fow of the video and 
positions of the UAV and the targets. 


Geo-referencing provides the synchronization between 3 dis- 
play modes: the 2D map, the 3D environment and the radar 
representation (graphical representation is discussed in the 
next subsection). It synchronizes in real-time the trajectory 
of the aircraft, the relative positions of possible targets and 
the actual video f ow on different types of maps. 


Payloads unit includes all payloads except camera. It in- 
cludes simulation of noise and bias on payloads. 


The Targets management unit allows the user to manage the 
fxed or mobile targets (coordinates, speed, trajectory, etc.). 


Graphical representation 


3D representation—A brief review of the signif cant features 
of 14 simulators (Flightgear, X-Plane, Simbad, Matlab, etc.) 
is presented in paper [11]. In order to evaluate the simulators 
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c Multiplay = in 
IP address PC1 
Port3 


State IP address PC2 
target Port 2 


Multiplay = out 


IP address PC1 
Port3 


Figure 2. Flightgear conf guration 


from the perspective of a robotics researcher, the authors 
propose four criteria: 


+ Physical Fidelity 

. Functional Fidelity 

« Ease of Development 

e Cost 

On the basis of this evaluation, we decided to work with Mat- 
lab/ Simulink and Flightgear. Matlab/Simulink is a numerical 
simulation environment which can communicate with other 
simulation environments that may provide better physics 
simulation and visualization. Also, it offers opportunity to 
work with the S-function language that allows programming 
separately the different modules and it is able to produce 
ordinary C programs. 


Flightgear has been chosen to simulate in real-time the virtual 
environment and the video fow. It is open-source so that the 
code can be adapted as needed. It runs on many different 
operating systems (Microsoft Windows and Linux) and it 
interfaces nicely with Matlab/Simulink. 


To simulate multiple objects (camera, targets, etc.) in the 
same Flightgear’s instance, servers or other computers are 
required. Hence it is possible to run multiple instances 
of Flightgear, the counterpart being the management of 
computer communications. In our case we have used two 
synchronized computers and we have added the Multiplay’s 
function to the script of Flightgear. As shown in Fig.2, we 
have defned the receiving PC (Multiplay = in) and sending 
PC (Multiplay = out) with corresponding port number. 


2D representation—The UAV trajectory is displayed on a 2D 
map (see Fig.3). Of course, connection and synchronization 
of the map with the Flightgear environment must be en- 
sured. For this purpose, we could have used the open-source 
software Atlas. Actually the choice between open-source 
or commercial software can be discussed. Atlas viewing is 
open-source application which can be connected directly to 
Flightgear. This application is a simple and free 2D map; it 
allows displaying aircraft’s current location on a moving map 
(121, [13]. 

However, Atlas is limited compared to our needs: 1. dynamic 
interaction with the map (def nition/redef nition of waypoints 
by clicking directly on the map); 2. display orientation of the 
camera’s line of sight; etc. 


To meet these requirements, we have chosen to develop our 
own 2D map application on Matlab/Simulink. For this we 
have designed a map by using aviation charts [14] to display 
the aircraft trajectory and waypoints. 
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Figure 3. 2D and 3D graphical representation 


On the other hand, we wanted to incorporate a radar view in 
our simulator, providing information about the relative posi- 
tions of possible targets. A schematic showing relationships 
between the set of display screens is depicted in Fig.3. 

The S-function “Bridge” ensures the communication between 
Matlab/Simulink and Flightgear. 


4. MAN MACHINE INTERFACE 


The set of simulator’s features are managed thanks to the 
Man Machine Interface (MMI) as shown in Fig.4. The user 
can therefore control and assess the state of the system. 
This interface deals with the information fow between the 
operator and UAV. The graphical interface proposes several 
windows which allow to set the parameters of simulation and 
to visualize the state of the system in real time. 


The MMI allows also initializing the simulator. This initial- 
ization phase is required before any operation. It permits 
to defne parameters that are specif c to a given simulation 
instance. The user selects a UAV type (f xed or rotary wing) 
among a data base. The aerodynamic model, the correspond- 
ing aerodynamic parameters and the model of payload are 
loaded accordingly. 


Finally the set of maps corresponding to the area where 
the scenario takes place, and the initial state of the aircraft 
are selected. Fig.4 shows a part of our MMI. The user 
can conf gure online simulation according to a predef ned 
scenario according to the following parameters: 


+ Waypoints table for the fight plan. The table can be 
completed in two distinct ways: entering the coordinates 
manually or directly clicking on the 2D map to automatically 
retrieve the coordinates of the points 

+ Navigation mode (navigation by waypoints, tracking, pat- 
tern of monitoring, etc.) 

e UAV control mode and camera control mode (manual or 
automatic) 

+ Geometrical patterns (circle radius, length and orientation 
of the racetrack, etc.) 

- Target management (position in case of a f xed target, speed 
and direction when a moving target) 

- Weather conditions: speed and direction of the wind. 


At this point, let us recall that one of our main purposes is to 
give the UAV a certain degree of autonomy (replanif cation). 
For instance when the initial mission is interrupted and when 
there is the need to determine a new path online, the user 
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can invoke certain replanif cation algorithms. Our algorithms 
are mostly based on Lyapunov-type method or on the Pon- 
tryagin Maximum Principle (PMP) [15]. We have chosen to 
minimize either the time or a length-curvature cost. Manual 
replanif cation can also be achieved by simply clicking on 
the map to create new waypoints. The choice between the 
different solutions is performed via MMI. 


At the level of the user the MMI displays several simulation 
pieces of information, such as: 


Flight data (position, heading, speed, etc.) 
3D video Flow 

UAV 2D trajectory 

Radar display 

Waypoints list 


5. AUTOMATIC CONTROL 


In this section, a description of the various control modules 
is proposed. The automation level is more or less signif cant 
depending on the involvement degree of the operator in the 
decision process. The control action can be performed by 
humans or computers. Thus, according to [16], [17], various 
levels of automation can be incorporated in a decision making 
system. 


The purpose of the control system at the lowest level is to 
generate control signals based upon position and attitude of 
the aircraft in order to satisfy the requirements of longitudinal 
and lateral maneuverability. 


The control module dedicated to the fxed wing model is 
depicted in Fig.5. This module is divided into four distinct 
hierarchical layers: fight path, navigation controller, autopi- 
lot and fight control surfaces controller. The top level of the 
control architecture manages the fight path. It encapsulates 
all the algorithms that generate a set of waypoints P, — 
(Ën; Va, Tn)n>0 with 


ên € R° : world coordinates of the point to be reached 
V, € 
T, € 


T : desired speed between P, and P,-: 


T : latency above the point (in the case of a helicopter) 


The navigation controller transforms the waypoints P,, into a 
desired velocity command VA, a desired altitude command h4 
and a desired heading command 4. 


The autopilot receives these commands and calculates head- 
ing, airspeed and altitude holders. The longitudinal and 
lateral dynamics of the carrier are considered decoupled so 
that the longitudinal and lateral autopilots are independent 


[6]. 


Using the measured altitude, speed w and angular velocity 
wo, the longitudinal autopilot gives the elevator def ection. 
The lateral autopilot depends on the desired heading, pitch 
angle ®, and angular velocity w.. The output is the heading 
command. 

A classical PID controller is enough for each autopilot. 


Finally, level zero consists of a direct access to the control 
surfaces of the UAV. Although this task is unrealistic we kept 
this possibility for internal tests (or for the fun). 

User can choose among three different modes for both the 
UAV and the camera turret: 
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Figure 4. The Man Machine Interface of the simulator 
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Figure 5. Level based decomposition of the controller 


Manual mode. A joystick allows to directly control the head- 
ing, speed and altitude commands. The azimuth and elevation 
commands of the camera can also be controlled. 
Semi-manual mode. In this case, at least one of the three 
components of the heading, speed and altitude command is 
managed by the system. 

Automatic mode. The UAV is fully autonomous and the sys- 
tem manages automatically the missions specif ed by the user. 


6. PATH PLANNING 


In this section, we present brief y our methodology for the 
trajectory planif cation/replanif cation for fxed wing UAV 
(Again, the rotary wing case is more classical and in fact 
easier, and we have developed a general methodology in a 
series of papers [18], [191, [20], [211, [221], [23]). 

In fxed wing case, we present two path planning approaches 
that can be used in replanif cation in order to join a pattern or 
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to track targets. The results will be stated without proof since 
they will be the subject of a forthcoming paper. The tangent 
vector feld guidance method and the Lyapunov method are 
used for UAV”Ss fight guidance. 

These approaches are limited by different constraints such 
as the geographical conf guration, the UAV’s dynamic char- 
acteristics, the UAV'Ss initial orientation and the traveling 
direction of the pattern. 

Finally, we propose to use an extended Kalman flter for 
target motion prediction. 


Tangent vector f eld guidance 


We present a dynamic path planning algorithm for a UAV 
based on the tangent vector feld guidance law that allows 
monitoring and tracking a ground target. This algorithm cal- 
culates a path by taking into account the physical constraints, 
the initial position of UAV and the position of target [241]. 


In Fig.6, UAV and target are represented in a topographical 
view. The UAV direction vector is denoted V(#). The goal 
is to join tangentially a circle with radius r around the target 
from any initial direction of the UAV. 


From the position of UAV we calculate two tangents to the 
target circle T'1 and T2. 1", 2 are the angles between 


the direction vector V(“) and T1, T2 respectively. In the 
case where the traveling direction along the f nal pattern is 
not def ned, the UAV should head the nearest tangent. An 
example is depicted in Fig.7. 

The new UAV direction angle 4 can be determined as: 


Va = min(y1, V2) 


The particular case when ||D|| < r is treated as follows: we 
simply allow the UAV to go straight ahead maintaining its 
initial direction. When the UAV is outside the circle, we apply 
again the algorithm. Notice that when the UAV starts far from 
the target, this strategy coincides exactly with the minimum 
time optimal policy [25]. 


target." 


Figure 6. From the position of UAV we calculate two 
tangents on the circle of radius r around the target 
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Figure 7.  Trajectory resulting from the tangent vector 
feld guidance starting at the point (15:0) with a direction 
of 0 radian and arriving tangent to a pattern, here the circle 
centering at (0,0) 


Lyapunoy vector f eld guidance 


In this section we consider the problem from the kinematic 
point of view only. In particular we consider that the aircraft 
commands are identically equal to angular velocities, 1.e 
(w1,wa,w3) = (ôl,0m,dn). We deal successively with the 
2D and 3D problem. The 2D problem is for instance the case 
of a HALE UAV at constant altitude. 


The Lyapunov method has been used by several authors, in 
several different ways [24], [26]. Here we present an original 
Lyapunov strategy, serving the advantage to be completely 
smooth, and very simple to apply. Frequently Lyapunov 
strategies presented in the litterature are discontinuous. 


First 2D problem—Our goal is to stabilize the UAV on a 
predef ned pattern (a circle) in the f ying plan above the target. 
From the kinematic point of view, the equations of the motion 
are just the standard Dubins equations (see [27] for a justif - 
cation), 1.e. 


à = cos 
ÿ =sinY (7) 
À =ws, 


with q = (x,y,v) € R? x S! being the state (where (x, y) € 
R? is the UAV's coordinates vector in the constant altitude 
plane, being the yaw angle), and w3 being the control 
variable. The yaw control w3 is constrained by wWazmin < 
w3 < W3 max» with W3 max > 0 and W3 max > W3 min: 


And we fx the target pattern to be a circle C: 
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Figure 8. LaSalle’s principle based planif cation results in 
2D 
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Figure 9. LaSalle’s principle based planif cation results in 
3D 


. centered at the origin with radius r = 1/43 max; 
- traveled counterclockwise. 


Let us observe that C corresponds to the maximal curvature 
that can be achieved by the system. 


We have the following result. Proof will be develop in a 
forthcoming paper. 


Theorem 6.1: Therexists an explicit feedback control func- 
tion u(-) which has the following properties: 


e u(-)is C®; 


. The pattern C is a globaly asymptotically stable attractor 
for the closed-loop system that result from applying u(-). 


Remark 6.2: In other Lyapunov type methods there are in 
general the 2 following drawbacks: 


1. The synthesis is not smooth 


2. In our method the pattern which is asymptotically reached 
is exactly the pattern C. In general it is not the case. For other 
Lyapunov technique f nal curvature is a bit smaller than the 
maximum curvature. 


Secondly 3D problem—Now we consider a MALE UAV 
fying at constant speed. Our goal does not change: we want 
to stabilize the UAV on a predef ned pattern above the target. 


Equations of motion are def ned on section 2. We considered 
that the UAV is controlled by three commands on angular 
velocity acting physically on three angles: roll, pitch and yaw 
def ned by: 


° wj, the roll control constrained by wi min < 1 < W1 max» 


with w1 max > 0 and W1 max >. W1 mins 
e wo, the pitch control constrained by wa min < W2 < W9 maxs 
with W9 max > 0 and W9 max > W9 min; 
° W3, the yaw control constrained by w3 min < W3 < W3 max, 
With W3 max > 0 and w3 max > 3 min: 


Here, we f x the pattern as the circle C3p: 
. centered at the origin with radius r — 


contained in a horizontal plane; 
. traveled counterclockwise. 


1/43 max and 


As for the 2D case we have the following result. 


Theorem 6.3: There exists three explicit control functions 
usp(r), à € {1,2,3}, for each UAV's controls, which have 
the following properties: 


. uip(-)is C®, i € {1,2,3}; 
. The pattern C3p is a globaly asymptotically stable attractor 
for the closed-loop system that result from applying the three 


control functions. 


Remark 6.4: Again one strong point is that the attractor C3 
is exactly reached. 


60 ———— UAV's path 
target's path 


50 100 150 


Figure 10. LaSalle’s principle based tracking results in 2D 
when the target velocity is close to the UAV”s velocity 
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Figure 11. LaSalle’s principle based tracking results in 2D 
when the target velocity is much smaller than UAV”Ss velocity 


Numerical results — We show numerical simulation using 
control functions of both Theorem 6.1 and 6.3. 
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First we present an example of a Dubin’s system simulation 
using the feedback control law of Theorem 6.1. Fig.8 shows 
a trajectory starting at the point (20,0) with a direction of 0 
radian and reaching circle C. 


Numerical result of Theorem 6.3 is showned on Fig.9. Here 
the UAV starts at the point (0,0,3) with an orientation of 0 
radian for the pitch, roll and yaw. The trajectory reach the 
circle C3p. 


Tracking of a moving target-—Algorithm shown in previous 
section gives some results when the target is f xed in location. 
Based on this algorithm we have developed another one 
algorithm to track a moving target. The idea is to readjust 
the pattern’s position according to the target’s position which 
can be known or estimated. 

To estimate the target position (when the target trajectory is 
unknown) we use an extended Kalman flter. It predicts the 
target velocity vector and estimates its future position. 


The simple case of a HALE UAV tracking a target with a 
priori known movement is shown on Fig.10 and 11. The 
UAV starts horizontaly at the position (10, —20,0) and the 
target at the origin. We present two cases depending on 
the difference between the target’s velocity and the aircraft’s 
velocity. Fig.10 shows result when target’s velocity is close 
to the aircraft’s velocity and Fig.11 shows a trajectory of a 
target having a lower speed. 

Fig.12 shows result of tracking a target following an hippo- 
drome pattern. The UAV starts horizontaly at the position 
(20; 30;10). In this example we use an extended Kalman 
f Iter to estimate the future position of the target. 


———— UAV's path 
10 target's path 
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Figure 12. LaSalle’s principle based tracking results in 3D 


7. SIMULATION RESULTS 


In this section we present simulation results of a fight sce- 
nario for a fxed wing UAV on our simulator. The assump- 
tions are: 


e the UAV has a suff cient reserve of fuel 

e the UAV moves at a f xed speed and constant altitude 

e the speed of the target is assumed to be f xed and less than 
the minimum speed of the UAV 

. wind speed is fxed at 5 m/s throughout the simulations. 


The Fig.13 depicts a fight scenario which implements several 
phases of the path planning and replanif cation. 


Figure 13. Diagram of complete scenario 


The scenario is as follows: 


A — B: (A) is the starting point. The UAV travels through 
the waypoints of a fight plan previously prepared. 

A1 — C: Before arriving at destination (B), the operator 
decides to divert the UAV and determines a new point of 
interest (C). In practice, the user clicks on the map or specif es 
manually coordinates of (C). The operator selects the replani- 
f cation algorithm to accomplish this task, knowing that there 
is a no-f y-zone on the way. 

C'— D: The operator decides to go to watch a stationary 
target located at (D). He chooses to turn around it by realizing 
a circle pattern. 

D — E: At (D) a moving target is identif ed, moving in the 
direction of (E), the operator decides to track it. 

E — A: The operator decides to end the mission and to 
return the UAV to its home base. 


Below, we present the simulation results by detailing each 
step. 


A — B: See Fig.14. The UAV starts its fight with an initial 
heading of 0 radians and following 3 waypoints spaced 2 
km at a speed of 80 knots. It passes through these points 
straightly. 


Figure 14. UAV trajectory between point À and C. UAV 
trajectory resulting of navigation by waypoint (A — A1) and 
Pontryagin Maximum Principle (A1 — C) 


A1 — C!: (Fig.14 again) At (A1) the UAV changes direction 
and moves to (C). To f nd an optimal trajectory and to avoid 
the obstacle on the way, we use a path planning algorithm 
based on the Pontryagin Maximum Principle (PMP). The 
PMP allows computing an optimal trajectory with a mini- 
mization of a cost function which is in our case the trade-off 
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between curve length and curvature. 


C'— D: In Fig.15 the purpose is to go monitoring a fxed 
target located at (D). For this we use the algorithm of tangent 
vector f eld guidance given in section 6 allowing the UAV to 
reach a pattern tangentially and we start describing a circle 
centered at (D) with a radius of 1km. 
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Figure 15. UAV trajectory resulting of the tangent vector 
feld guidance starting at point C and arriving tangentially at 
the pattern centered at point D. 


Figure 16. Camera view of the f xed target. 


D — E: Fig.17 shows the trajectory of the UAV tracking a 
moving target which moves straight at a constant speed 45 
knots and less the minimum speed of the UAV which moves 
at 70 knots (its minimum speed). To achieve this mission we 
use the tracking algorithm based on the Lyapunov vector f eld 
guidance algorithm presented in Section 6. 
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Figure 17. UAV trajectory resulting of the Lyapunov vector 
feld guidance starting at point D and tracking a target moving 
in a straight line (red line) at a constant speed. 
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Figure 18. Path back to the base following the waypoints 
def ned online 


E — À: In Fig.18, the UAV returns back in navigation mode 
and follows waypoints set to return to the base. 


These simulations allowed us to test and validate our simula- 
tor (model, control module, MMI, synchronization between 
different mode of view, etc.) through a reasonable scenario 
and to test our replanif cation algorithms. 


8. CONCLUSION 


In this paper we presented a ground control station simulator 
for fxed and rotary wing. The simulator’s architecture is 
based upon modules programmed independently, allowing 
modularity and easy maneuverability. As a consequence the 
simulator 1s able to simulate any UAV model. 


Besides the well known advantages of the modularity, the 
operator can select simulation parameters and visualize the 
behavior of UAV and target in a virtual environment, thanks 
to MMI interface. 

The simulator allowed us to simulate and validate all the 
operational environment of the station. Furthermore, the 
simulator could be used as a portable demonstrator, or for 
the purpose of operators training. 


Another part of this paper concerns the path planif ca- 
tion/replanif cation algorithms for fxed wing UAVS. Actu- 
ally, our simulated station is assumed to provide a certain 
degree of autonomy to the UAV. 

We developed certain original path planning algorithms in 2D 
and 3D. 


Finally, we presented simulation results for scenarios of 
missions. To conclude this work we display in Fig.19 a view 
of the ground control station simulator. 


9. FUTURE WORKS 


One of the goals of this work is to develop a ground control 
station simulator and to def ne eff cient algorithms for tra- 
jectory planif cation/replanif cation. In previous subsections, 
we have presented two path planning algorithms and an 
Extended Kalman Filter has been chosen for the target motion 
estimation. 


One of the future work directions concerns the improvement 
of estimated location of the target. In this context, the idea is 
to improve on target tracking with road network information. 
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Many studies [24], [28], [29] deal with this problem and 
different approaches are proposed. For the target motion 
prediction the authors use particle fIters and a motion target 
model can be used depending on the target characteristics. 


Another idea is to develop a decision support system in order 
to aid the operator in a replanif cation context of mission. 
This system should allow the operator to take the best deci- 
sion under a set of rules/constraints. To propose a new path, 
this system should take into account safety information like 
fuel rate, weather conditions, risk zones, time objectives and 
environment parameters. Several approaches are possible for 
this purpose, such as fuzzy logic, bayesian networks, neural 
networks or multicriteria methods [30], [31], [32], [33]. 
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ANNEXE B. ARTICLE DE CONFÉRENCE CONCERNANT L’APPLICATION DES 
MÉTHODES DE PLANIFICATION SUR LE SIMULATEUR [?] : “ATH PLANNING AND 
GROUND CONTROL STATION SIMULATOR FOR UAV” 
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Title : Path planning of unmanned combat aircraft vehicles 


Abstract : This thesis is about path planning for HALE or MALE UAVSs (Unmanned Aircraft Vehicles), 
possibly under mission constraints. As such, the study is performed at the kinematic level : HALE UAVs are 
represented as Dubins systems, and à model for MALE UAVSs is constructed by studying their kinematic frame. 

In the first part, we tackle the path planning problem for a UAV that must join a target (a point or a pattern), 
starting from any position. The point to point path planning problem is addressed as an optimal control problem. 
Regarding the point to pattern path planning problem, two different methods are proposed. The former consists in 
solving the minimum time synthesis for the Dubins system, in order to obtain à basis for a HALE UAVSs planning 
algorithm. The latter method relies on the LaSalle principle; it permits to stabilize a HALE or MALE UAV to a 
pattern. 

In addition, extensions of the previously developed algorithms to cluttered environnement are provided. This 
extension is achieved thanks to a space discretization using Voronoi diagrams and a discrete planning method. 

Finally, the mission constraints are dealt with as à coupling problem between the UAV and its sensors. The 
proposed algorithm is presented in the form of a constrained quadratic problem. 

In the second part of this thesis, we want to refine the planning algorithm to get a result closer to trajectories 
of pilots. In order to do that, we solve an inverse optimal control problem where the cost to find is computed 
from the experience of pilots. Theoretical results are presented and applied to the particular case of the Dubins 
system. 


Keywords : Optimal control, Lyapunov method, LaSalle principle, inverse optimal control problem, anthropo- 
morphic control, drone, identification, modeling, path planning. 


Titre : Planification de trajectoire pour drones de combat 


Résumé : L'objectif principal de ce travail est l’étude de la planification de trajectoires pour des drones de 
type HALE ou MALE. Les modèles cinématiques de ces drones sont étudiés. Les drones HALE sont modélisés 
par le système de Dubins. Pour les drones MALE, le modèle est construit en étudiant le repère cinématique du 
drone. 

Nous considérons les problèmes de planification de trajectoires point-point et point-pattern. Il s’agit, à partir 
de la position courante du drone, de rejoindre un point ou une figure prédéfinie dans l’espace. La planification 
point-point est abordée sous forme d’un problème de contrôle optimal. Deux méthodes sont proposées pour 
résoudre le problème point-pattern. D'abord nous présentons la synthèse en temps minimal pour le système de 
Dubins. Ensuite, nous développons une méthode basée sur le principe de LaSalle. La première méthode est utilisée 
au sein d’un algorithme de planification pour des drones HALE. La deuxième permet de stabiliser les deux types 
de drones considérés vers un pattern. 

Nous proposons une extension des algorithmes de planification développés, basée sur une discrétisation de 
l’espace grâce aux graphes de Voronoï et une méthode de planification discrète, pour construire des trajectoires 
dans des milieux encombrés. 

Nous étudions également le problème de couplage drone/capteur. Il s’agit de calculer une trajectoire permet- 
tant de satisfaire les objectifs du drone et de son capteur (une caméra). L’algorithme proposé est construit à 
partir de la résolution d’un problème quadratique sous contraintes. 

Dans une seconde partie, nous analysons un problème de contrôle optimal inverse. Celui-ci permet d’amélio- 
rer les résultats des méthodes de planification en s'inspirant du comportement des pilotes. Après avoir posé le 
problème, les résultats théoriques sont exposés et le cas particulier du système de Dubins est étudié en pratique. 


Mots-clefs : Contrôle optimal, méthode de Lyapunov, principe de LaSalle, problème de contrôle optimal inverse, 
commande anthropomorphique, drone, identification, modélisation, planification de trajectoires. 


